Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 : /*
11 : Private DS routines
12 : */
13 :
14 : #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/
15 : #include <slepcblaslapack.h>
16 :
17 20188 : PetscErrorCode DSAllocateMat_Private(DS ds,DSMatType m)
18 : {
19 20188 : PetscInt n,d,rows=0,cols;
20 20188 : PetscBool ispep,isnep;
21 :
22 20188 : PetscFunctionBegin;
23 20188 : n = ds->ld;
24 20188 : PetscCall(PetscObjectTypeCompare((PetscObject)ds,DSPEP,&ispep));
25 20188 : if (ispep) {
26 1136 : if (m==DS_MAT_A || m==DS_MAT_B || m==DS_MAT_W || m==DS_MAT_U || m==DS_MAT_X || m==DS_MAT_Y) {
27 1094 : PetscCall(DSPEPGetDegree(ds,&d));
28 1094 : n = d*ds->ld;
29 : }
30 : } else {
31 19052 : PetscCall(PetscObjectTypeCompare((PetscObject)ds,DSNEP,&isnep));
32 19052 : if (isnep) {
33 635 : if (m==DS_MAT_Q || m==DS_MAT_Z || m==DS_MAT_U || m==DS_MAT_V || m==DS_MAT_X || m==DS_MAT_Y) {
34 20 : PetscCall(DSNEPGetMinimality(ds,&d));
35 20 : n = d*ds->ld;
36 : }
37 : }
38 : }
39 20188 : cols = n;
40 :
41 20188 : switch (m) {
42 : case DS_MAT_T:
43 : cols = PetscDefined(USE_COMPLEX)? 2: 3;
44 : rows = n;
45 : break;
46 151 : case DS_MAT_D:
47 151 : cols = 1;
48 151 : rows = n;
49 151 : break;
50 4074 : case DS_MAT_X:
51 4074 : rows = ds->ld;
52 4074 : break;
53 81 : case DS_MAT_Y:
54 81 : rows = ds->ld;
55 81 : break;
56 15349 : default:
57 15349 : rows = n;
58 : }
59 20188 : if (ds->omat[m]) PetscCall(MatZeroEntries(ds->omat[m]));
60 : else {
61 4354 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,rows,cols,NULL,&ds->omat[m]));
62 : }
63 20188 : PetscFunctionReturn(PETSC_SUCCESS);
64 : }
65 :
66 91776 : PetscErrorCode DSAllocateWork_Private(DS ds,PetscInt s,PetscInt r,PetscInt i)
67 : {
68 91776 : PetscFunctionBegin;
69 91776 : if (s>ds->lwork) {
70 1508 : PetscCall(PetscFree(ds->work));
71 1508 : PetscCall(PetscMalloc1(s,&ds->work));
72 1508 : ds->lwork = s;
73 : }
74 91776 : if (r>ds->lrwork) {
75 505 : PetscCall(PetscFree(ds->rwork));
76 505 : PetscCall(PetscMalloc1(r,&ds->rwork));
77 505 : ds->lrwork = r;
78 : }
79 91776 : if (i>ds->liwork) {
80 686 : PetscCall(PetscFree(ds->iwork));
81 686 : PetscCall(PetscMalloc1(i,&ds->iwork));
82 686 : ds->liwork = i;
83 : }
84 91776 : PetscFunctionReturn(PETSC_SUCCESS);
85 : }
86 :
87 : /*@C
88 : DSViewMat - Prints one of the internal DS matrices.
89 :
90 : Collective
91 :
92 : Input Parameters:
93 : + ds - the direct solver context
94 : . viewer - visualization context
95 : - m - matrix to display
96 :
97 : Note:
98 : Works only for ascii viewers. Set the viewer in Matlab format if
99 : want to paste into Matlab.
100 :
101 : Level: developer
102 :
103 : .seealso: DSView()
104 : @*/
105 28 : PetscErrorCode DSViewMat(DS ds,PetscViewer viewer,DSMatType m)
106 : {
107 28 : PetscInt i,j,rows,cols;
108 28 : const PetscScalar *M=NULL,*v;
109 28 : PetscViewerFormat format;
110 : #if defined(PETSC_USE_COMPLEX)
111 : PetscBool allreal = PETSC_TRUE;
112 : const PetscReal *vr;
113 : #endif
114 :
115 28 : PetscFunctionBegin;
116 28 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
117 112 : PetscValidLogicalCollectiveEnum(ds,m,3);
118 28 : DSCheckValidMat(ds,m,3);
119 28 : if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer));
120 28 : PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
121 28 : PetscCheckSameComm(ds,1,viewer,2);
122 28 : PetscCall(PetscViewerGetFormat(viewer,&format));
123 28 : if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
124 28 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
125 28 : PetscCall(DSMatGetSize(ds,m,&rows,&cols));
126 28 : PetscCall(MatDenseGetArrayRead(ds->omat[m],&M));
127 : #if defined(PETSC_USE_COMPLEX)
128 : /* determine if matrix has all real values */
129 : if (m!=DS_MAT_T && m!=DS_MAT_D) {
130 : /* determine if matrix has all real values */
131 : v = M;
132 : for (i=0;i<rows;i++)
133 : for (j=0;j<cols;j++)
134 : if (PetscImaginaryPart(v[i+j*ds->ld])) { allreal = PETSC_FALSE; break; }
135 : }
136 : if (m==DS_MAT_T) cols=3;
137 : #endif
138 28 : if (format == PETSC_VIEWER_ASCII_MATLAB) {
139 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols));
140 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s = [\n",DSMatName[m]));
141 28 : } else PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[m]));
142 :
143 388 : for (i=0;i<rows;i++) {
144 360 : v = M+i;
145 : #if defined(PETSC_USE_COMPLEX)
146 : vr = (const PetscReal*)M+i; /* handle compact storage, 2nd column is in imaginary part of 1st column */
147 : #endif
148 3840 : for (j=0;j<cols;j++) {
149 : #if defined(PETSC_USE_COMPLEX)
150 : if (m!=DS_MAT_T && m!=DS_MAT_D) {
151 : if (allreal) PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v)));
152 : else PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v)));
153 : } else PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*vr));
154 : #else
155 3480 : PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v));
156 : #endif
157 3480 : v += ds->ld;
158 : #if defined(PETSC_USE_COMPLEX)
159 : if (m==DS_MAT_T) vr += ds->ld;
160 : #endif
161 : }
162 360 : PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
163 : }
164 28 : PetscCall(MatDenseRestoreArrayRead(ds->omat[m],&M));
165 :
166 28 : if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
167 28 : PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
168 28 : PetscCall(PetscViewerFlush(viewer));
169 28 : PetscFunctionReturn(PETSC_SUCCESS);
170 : }
171 :
172 20065 : PetscErrorCode DSSortEigenvalues_Private(DS ds,PetscScalar *wr,PetscScalar *wi,PetscInt *perm,PetscBool isghiep)
173 : {
174 20065 : PetscScalar re,im,wi0;
175 20065 : PetscInt n,i,j,result,tmp1,tmp2=0,d=1;
176 :
177 20065 : PetscFunctionBegin;
178 20065 : n = ds->t; /* sort only first t pairs if truncated */
179 : /* insertion sort */
180 20065 : i=ds->l+1;
181 : #if !defined(PETSC_USE_COMPLEX)
182 20065 : if (wi && wi[perm[i-1]]!=0.0) i++; /* initial value is complex */
183 : #else
184 : if (isghiep && PetscImaginaryPart(wr[perm[i-1]])!=0.0) i++;
185 : #endif
186 278972 : for (;i<n;i+=d) {
187 258907 : re = wr[perm[i]];
188 258907 : if (wi) im = wi[perm[i]];
189 : else im = 0.0;
190 258907 : tmp1 = perm[i];
191 : #if !defined(PETSC_USE_COMPLEX)
192 258907 : if (im!=0.0) { d = 2; tmp2 = perm[i+1]; }
193 : else d = 1;
194 : #else
195 : if (isghiep && PetscImaginaryPart(re)!=0.0) { d = 2; tmp2 = perm[i+1]; }
196 : else d = 1;
197 : #endif
198 258907 : j = i-1;
199 258907 : if (wi) wi0 = wi[perm[j]];
200 : else wi0 = 0.0;
201 258907 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result));
202 1252181 : while (result<0 && j>=ds->l) {
203 1170111 : perm[j+d] = perm[j];
204 1170111 : j--;
205 : #if !defined(PETSC_USE_COMPLEX)
206 1170111 : if (wi && wi[perm[j+1]]!=0)
207 : #else
208 : if (isghiep && PetscImaginaryPart(wr[perm[j+1]])!=0)
209 : #endif
210 47062 : { perm[j+d] = perm[j]; j--; }
211 :
212 1170111 : if (j>=ds->l) {
213 1088041 : if (wi) wi0 = wi[perm[j]];
214 : else wi0 = 0.0;
215 2517059 : PetscCall(SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result));
216 : }
217 : }
218 258907 : perm[j+1] = tmp1;
219 258907 : if (d==2) perm[j+2] = tmp2;
220 : }
221 20065 : PetscFunctionReturn(PETSC_SUCCESS);
222 : }
223 :
224 5867 : PetscErrorCode DSSortEigenvaluesReal_Private(DS ds,PetscReal *eig,PetscInt *perm)
225 : {
226 5867 : PetscScalar re;
227 5867 : PetscInt i,j,result,tmp,l,n;
228 :
229 5867 : PetscFunctionBegin;
230 5867 : n = ds->t; /* sort only first t pairs if truncated */
231 5867 : l = ds->l;
232 : /* insertion sort */
233 90130 : for (i=l+1;i<n;i++) {
234 84263 : re = eig[perm[i]];
235 84263 : j = i-1;
236 84263 : PetscCall(SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result));
237 513909 : while (result<0 && j>=l) {
238 487307 : tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
239 1058877 : if (j>=l) PetscCall(SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result));
240 : }
241 : }
242 5867 : PetscFunctionReturn(PETSC_SUCCESS);
243 : }
244 :
245 : /*
246 : Permute comumns [istart..iend-1] of [mat] according to perm. Columns have length n
247 : */
248 22711 : PetscErrorCode DSPermuteColumns_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat,PetscInt *perm)
249 : {
250 22711 : PetscInt i,j,k,p,ld;
251 22711 : PetscScalar *Q,rtmp;
252 :
253 22711 : PetscFunctionBegin;
254 22711 : ld = ds->ld;
255 22711 : PetscCall(MatDenseGetArray(ds->omat[mat],&Q));
256 350907 : for (i=istart;i<iend;i++) {
257 328196 : p = perm[i];
258 328196 : if (p != i) {
259 152240 : j = i + 1;
260 1257133 : while (perm[j] != i) j++;
261 152240 : perm[j] = p; perm[i] = i;
262 : /* swap columns i and j */
263 4501195 : for (k=0;k<n;k++) {
264 4348955 : rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
265 : }
266 : }
267 : }
268 22711 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&Q));
269 22711 : PetscFunctionReturn(PETSC_SUCCESS);
270 : }
271 :
272 : /*
273 : The same as DSPermuteColumns_Private but for two matrices [mat1] and [mat2]
274 : */
275 267 : PetscErrorCode DSPermuteColumnsTwo_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat1,DSMatType mat2,PetscInt *perm)
276 : {
277 267 : PetscInt i,j,k,p,ld;
278 267 : PetscScalar *Q,*Z,rtmp,rtmp2;
279 :
280 267 : PetscFunctionBegin;
281 267 : ld = ds->ld;
282 267 : PetscCall(MatDenseGetArray(ds->omat[mat1],&Q));
283 267 : PetscCall(MatDenseGetArray(ds->omat[mat2],&Z));
284 5813 : for (i=istart;i<iend;i++) {
285 5546 : p = perm[i];
286 5546 : if (p != i) {
287 4021 : j = i + 1;
288 38032 : while (perm[j] != i) j++;
289 4021 : perm[j] = p; perm[i] = i;
290 : /* swap columns i and j */
291 56518 : for (k=0;k<n;k++) {
292 52497 : rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
293 52497 : rtmp2 = Z[k+p*ld]; Z[k+p*ld] = Z[k+i*ld]; Z[k+i*ld] = rtmp2;
294 : }
295 : }
296 : }
297 267 : PetscCall(MatDenseRestoreArray(ds->omat[mat1],&Q));
298 267 : PetscCall(MatDenseRestoreArray(ds->omat[mat2],&Z));
299 267 : PetscFunctionReturn(PETSC_SUCCESS);
300 : }
301 :
302 : /*
303 : Permute rows [istart..iend-1] of [mat] according to perm. Rows have length m
304 : */
305 0 : PetscErrorCode DSPermuteRows_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt m,DSMatType mat,PetscInt *perm)
306 : {
307 0 : PetscInt i,j,k,p,ld;
308 0 : PetscScalar *Q,rtmp;
309 :
310 0 : PetscFunctionBegin;
311 0 : ld = ds->ld;
312 0 : PetscCall(MatDenseGetArray(ds->omat[mat],&Q));
313 0 : for (i=istart;i<iend;i++) {
314 0 : p = perm[i];
315 0 : if (p != i) {
316 0 : j = i + 1;
317 0 : while (perm[j] != i) j++;
318 0 : perm[j] = p; perm[i] = i;
319 : /* swap rows i and j */
320 0 : for (k=0;k<m;k++) {
321 0 : rtmp = Q[p+k*ld]; Q[p+k*ld] = Q[i+k*ld]; Q[i+k*ld] = rtmp;
322 : }
323 : }
324 : }
325 0 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&Q));
326 0 : PetscFunctionReturn(PETSC_SUCCESS);
327 : }
328 :
329 : /*
330 : Permute columns [istart..iend-1] of [mat1] and [mat2] according to perm.
331 : Columns of [mat1] have length n, columns of [mat2] have length m
332 : */
333 2015 : PetscErrorCode DSPermuteBoth_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,PetscInt m,DSMatType mat1,DSMatType mat2,PetscInt *perm)
334 : {
335 2015 : PetscInt i,j,k,p,ld;
336 2015 : PetscScalar *U,*V,rtmp;
337 :
338 2015 : PetscFunctionBegin;
339 2015 : ld = ds->ld;
340 2015 : PetscCall(MatDenseGetArray(ds->omat[mat1],&U));
341 2015 : PetscCall(MatDenseGetArray(ds->omat[mat2],&V));
342 28061 : for (i=istart;i<iend;i++) {
343 26046 : p = perm[i];
344 26046 : if (p != i) {
345 5200 : j = i + 1;
346 75616 : while (perm[j] != i) j++;
347 5200 : perm[j] = p; perm[i] = i;
348 : /* swap columns i and j of U */
349 478939 : for (k=0;k<n;k++) {
350 473739 : rtmp = U[k+p*ld]; U[k+p*ld] = U[k+i*ld]; U[k+i*ld] = rtmp;
351 : }
352 : /* swap columns i and j of V */
353 277120 : for (k=0;k<m;k++) {
354 271920 : rtmp = V[k+p*ld]; V[k+p*ld] = V[k+i*ld]; V[k+i*ld] = rtmp;
355 : }
356 : }
357 : }
358 2015 : PetscCall(MatDenseRestoreArray(ds->omat[mat1],&U));
359 2015 : PetscCall(MatDenseRestoreArray(ds->omat[mat2],&V));
360 2015 : PetscFunctionReturn(PETSC_SUCCESS);
361 : }
362 :
363 : /*@
364 : DSSetIdentity - Copy the identity (a diagonal matrix with ones) on the
365 : active part of a matrix.
366 :
367 : Logically Collective
368 :
369 : Input Parameters:
370 : + ds - the direct solver context
371 : - mat - the matrix to modify
372 :
373 : Level: intermediate
374 :
375 : .seealso: DSGetMat()
376 : @*/
377 9610 : PetscErrorCode DSSetIdentity(DS ds,DSMatType mat)
378 : {
379 9610 : PetscScalar *A;
380 9610 : PetscInt i,ld,n,l;
381 :
382 9610 : PetscFunctionBegin;
383 9610 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
384 38440 : PetscValidLogicalCollectiveEnum(ds,mat,2);
385 9610 : DSCheckValidMat(ds,mat,2);
386 :
387 9610 : PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
388 9610 : PetscCall(DSGetLeadingDimension(ds,&ld));
389 9610 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
390 9610 : PetscCall(MatDenseGetArray(ds->omat[mat],&A));
391 9610 : PetscCall(PetscArrayzero(A+l*ld,ld*(n-l)));
392 137591 : for (i=l;i<n;i++) A[i+i*ld] = 1.0;
393 9610 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&A));
394 9610 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
395 9610 : PetscFunctionReturn(PETSC_SUCCESS);
396 : }
397 :
398 : /*@C
399 : DSOrthogonalize - Orthogonalize the columns of a matrix.
400 :
401 : Logically Collective
402 :
403 : Input Parameters:
404 : + ds - the direct solver context
405 : . mat - a matrix
406 : - cols - number of columns to orthogonalize (starting from column zero)
407 :
408 : Output Parameter:
409 : . lindcols - (optional) number of linearly independent columns of the matrix
410 :
411 : Level: developer
412 :
413 : .seealso: DSPseudoOrthogonalize()
414 : @*/
415 1606 : PetscErrorCode DSOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscInt *lindcols)
416 : {
417 1606 : PetscInt n,l,ld;
418 1606 : PetscBLASInt ld_,rA,cA,info,ltau,lw;
419 1606 : PetscScalar *A,*tau,*w,saux,dummy;
420 :
421 1606 : PetscFunctionBegin;
422 1606 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
423 1606 : DSCheckAlloc(ds,1);
424 6424 : PetscValidLogicalCollectiveEnum(ds,mat,2);
425 1606 : DSCheckValidMat(ds,mat,2);
426 6424 : PetscValidLogicalCollectiveInt(ds,cols,3);
427 :
428 1606 : PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
429 1606 : PetscCall(DSGetLeadingDimension(ds,&ld));
430 1606 : n = n - l;
431 1606 : PetscCheck(cols<=n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
432 1606 : if (n == 0 || cols == 0) PetscFunctionReturn(PETSC_SUCCESS);
433 :
434 1606 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
435 1606 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
436 1606 : PetscCall(MatDenseGetArray(ds->omat[mat],&A));
437 1606 : PetscCall(PetscBLASIntCast(PetscMin(cols,n),<au));
438 1606 : PetscCall(PetscBLASIntCast(ld,&ld_));
439 1606 : PetscCall(PetscBLASIntCast(n,&rA));
440 1606 : PetscCall(PetscBLASIntCast(cols,&cA));
441 1606 : lw = -1;
442 1606 : PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,A,&ld_,&dummy,&saux,&lw,&info));
443 1606 : SlepcCheckLapackInfo("geqrf",info);
444 1606 : lw = (PetscBLASInt)PetscRealPart(saux);
445 1606 : PetscCall(DSAllocateWork_Private(ds,lw+ltau,0,0));
446 1606 : tau = ds->work;
447 1606 : w = &tau[ltau];
448 1606 : PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,&A[ld*l+l],&ld_,tau,w,&lw,&info));
449 1606 : SlepcCheckLapackInfo("geqrf",info);
450 1606 : PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&rA,<au,<au,&A[ld*l+l],&ld_,tau,w,&lw,&info));
451 1606 : SlepcCheckLapackInfo("orgqr",info);
452 1606 : if (lindcols) *lindcols = ltau;
453 1606 : PetscCall(PetscFPTrapPop());
454 1606 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&A));
455 1606 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
456 1606 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
457 1606 : PetscFunctionReturn(PETSC_SUCCESS);
458 : }
459 :
460 : /*
461 : Compute C <- a*A*B + b*C, where
462 : ldC, the leading dimension of C,
463 : ldA, the leading dimension of A,
464 : rA, cA, rows and columns of A,
465 : At, if true use the transpose of A instead,
466 : ldB, the leading dimension of B,
467 : rB, cB, rows and columns of B,
468 : Bt, if true use the transpose of B instead
469 : */
470 55 : static PetscErrorCode SlepcMatDenseMult(PetscScalar *C,PetscInt _ldC,PetscScalar b,PetscScalar a,const PetscScalar *A,PetscInt _ldA,PetscInt rA,PetscInt cA,PetscBool At,const PetscScalar *B,PetscInt _ldB,PetscInt rB,PetscInt cB,PetscBool Bt)
471 : {
472 55 : PetscInt tmp;
473 55 : PetscBLASInt m, n, k, ldA = _ldA, ldB = _ldB, ldC = _ldC;
474 55 : const char *N = "N", *T = "C", *qA = N, *qB = N;
475 :
476 55 : PetscFunctionBegin;
477 55 : if ((rA == 0) || (cB == 0)) PetscFunctionReturn(PETSC_SUCCESS);
478 55 : PetscAssertPointer(C,1);
479 55 : PetscAssertPointer(A,5);
480 55 : PetscAssertPointer(B,10);
481 :
482 : /* Transpose if needed */
483 55 : if (At) tmp = rA, rA = cA, cA = tmp, qA = T;
484 55 : if (Bt) tmp = rB, rB = cB, cB = tmp, qB = T;
485 :
486 : /* Check size */
487 55 : PetscCheck(cA==rB,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Matrix dimensions do not match");
488 :
489 : /* Do stub */
490 55 : if ((rA == 1) && (cA == 1) && (cB == 1)) {
491 0 : if (!At && !Bt) *C = *A * *B;
492 0 : else if (At && !Bt) *C = PetscConj(*A) * *B;
493 0 : else if (!At && Bt) *C = *A * PetscConj(*B);
494 0 : else *C = PetscConj(*A) * PetscConj(*B);
495 0 : m = n = k = 1;
496 : } else {
497 55 : m = rA; n = cB; k = cA;
498 55 : PetscCallBLAS("BLASgemm",BLASgemm_(qA,qB,&m,&n,&k,&a,(PetscScalar*)A,&ldA,(PetscScalar*)B,&ldB,&b,C,&ldC));
499 : }
500 :
501 55 : PetscCall(PetscLogFlops(2.0*m*n*k));
502 55 : PetscFunctionReturn(PETSC_SUCCESS);
503 : }
504 :
505 : /*@C
506 : DSPseudoOrthogonalize - Orthogonalize the columns of a matrix with Modified
507 : Gram-Schmidt in an indefinite inner product space defined by a signature.
508 :
509 : Logically Collective
510 :
511 : Input Parameters:
512 : + ds - the direct solver context
513 : . mat - the matrix
514 : . cols - number of columns to orthogonalize (starting from column zero)
515 : - s - the signature that defines the inner product
516 :
517 : Output Parameters:
518 : + lindcols - (optional) linearly independent columns of the matrix
519 : - ns - (optional) the new signature of the vectors
520 :
521 : Note:
522 : After the call the matrix satisfies A'*s*A = ns.
523 :
524 : Level: developer
525 :
526 : .seealso: DSOrthogonalize()
527 : @*/
528 1 : PetscErrorCode DSPseudoOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscReal *s,PetscInt *lindcols,PetscReal *ns)
529 : {
530 1 : PetscInt i,j,k,l,n,ld;
531 1 : PetscBLASInt info,one=1,zero=0,rA_,ld_;
532 1 : PetscScalar *A,*A_,*m,*h,nr0;
533 1 : PetscReal nr_o,nr,nr_abs,*ns_,done=1.0;
534 :
535 1 : PetscFunctionBegin;
536 1 : PetscValidHeaderSpecific(ds,DS_CLASSID,1);
537 1 : DSCheckAlloc(ds,1);
538 4 : PetscValidLogicalCollectiveEnum(ds,mat,2);
539 1 : DSCheckValidMat(ds,mat,2);
540 4 : PetscValidLogicalCollectiveInt(ds,cols,3);
541 1 : PetscAssertPointer(s,4);
542 1 : PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
543 1 : PetscCall(DSGetLeadingDimension(ds,&ld));
544 1 : n = n - l;
545 1 : PetscCheck(cols<=n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
546 1 : if (n == 0 || cols == 0) PetscFunctionReturn(PETSC_SUCCESS);
547 1 : PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
548 1 : PetscCall(PetscBLASIntCast(n,&rA_));
549 1 : PetscCall(PetscBLASIntCast(ld,&ld_));
550 1 : PetscCall(MatDenseGetArray(ds->omat[mat],&A_));
551 1 : A = &A_[ld*l+l];
552 2 : PetscCall(DSAllocateWork_Private(ds,n+cols,ns?0:cols,0));
553 1 : m = ds->work;
554 1 : h = &m[n];
555 1 : ns_ = ns ? ns : ds->rwork;
556 11 : for (i=0; i<cols; i++) {
557 : /* m <- diag(s)*A[i] */
558 110 : for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
559 : /* nr_o <- mynorm(A[i]'*m), mynorm(x) = sign(x)*sqrt(|x|) */
560 10 : PetscCall(SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
561 10 : nr = nr_o = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
562 16 : for (j=0; j<3 && i>0; j++) {
563 : /* h <- A[0:i-1]'*m */
564 15 : PetscCall(SlepcMatDenseMult(h,i,0.0,1.0,A,ld,n,i,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
565 : /* h <- diag(ns)*h */
566 99 : for (k=0; k<i; k++) h[k] *= ns_[k];
567 : /* A[i] <- A[i] - A[0:i-1]*h */
568 15 : PetscCall(SlepcMatDenseMult(&A[ld*i],ld,1.0,-1.0,A,ld,n,i,PETSC_FALSE,h,i,i,1,PETSC_FALSE));
569 : /* m <- diag(s)*A[i] */
570 165 : for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
571 : /* nr_o <- mynorm(A[i]'*m) */
572 15 : PetscCall(SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
573 15 : nr = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
574 15 : PetscCheck(PetscAbs(nr)>PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linear dependency detected");
575 15 : if (PetscAbs(nr) > 0.7*PetscAbs(nr_o)) break;
576 6 : nr_o = nr;
577 : }
578 10 : ns_[i] = PetscSign(nr);
579 : /* A[i] <- A[i]/|nr| */
580 10 : nr_abs = PetscAbs(nr);
581 10 : PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&nr_abs,&done,&rA_,&one,A+i*ld,&ld_,&info));
582 10 : SlepcCheckLapackInfo("lascl",info);
583 : }
584 1 : PetscCall(MatDenseRestoreArray(ds->omat[mat],&A_));
585 1 : PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
586 1 : PetscCall(PetscObjectStateIncrease((PetscObject)ds));
587 1 : if (lindcols) *lindcols = cols;
588 1 : PetscFunctionReturn(PETSC_SUCCESS);
589 : }
|
