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