Actual source code: dspriv.c
slepc-3.22.1 2024-10-28
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 DSAllocateWork_Private(DS ds,PetscInt s,PetscInt r,PetscInt i)
67: {
68: PetscFunctionBegin;
69: if (s>ds->lwork) {
70: PetscCall(PetscFree(ds->work));
71: PetscCall(PetscMalloc1(s,&ds->work));
72: ds->lwork = s;
73: }
74: if (r>ds->lrwork) {
75: PetscCall(PetscFree(ds->rwork));
76: PetscCall(PetscMalloc1(r,&ds->rwork));
77: ds->lrwork = r;
78: }
79: if (i>ds->liwork) {
80: PetscCall(PetscFree(ds->iwork));
81: PetscCall(PetscMalloc1(i,&ds->iwork));
82: ds->liwork = i;
83: }
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: /*@
88: DSViewMat - Prints one of the internal DS matrices.
90: Collective
92: Input Parameters:
93: + ds - the direct solver context
94: . viewer - visualization context
95: - m - matrix to display
97: Note:
98: Works only for ascii viewers. Set the viewer in Matlab format if
99: want to paste into Matlab.
101: Level: developer
103: .seealso: DSView()
104: @*/
105: PetscErrorCode DSViewMat(DS ds,PetscViewer viewer,DSMatType m)
106: {
107: PetscInt i,j,rows,cols;
108: const PetscScalar *M=NULL,*v;
109: PetscViewerFormat format;
110: #if defined(PETSC_USE_COMPLEX)
111: PetscBool allreal = PETSC_TRUE;
112: const PetscReal *vr;
113: #endif
115: PetscFunctionBegin;
118: DSCheckValidMat(ds,m,3);
119: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer));
121: PetscCheckSameComm(ds,1,viewer,2);
122: PetscCall(PetscViewerGetFormat(viewer,&format));
123: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
124: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
125: PetscCall(DSMatGetSize(ds,m,&rows,&cols));
126: 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: if (format == PETSC_VIEWER_ASCII_MATLAB) {
139: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols));
140: PetscCall(PetscViewerASCIIPrintf(viewer,"%s = [\n",DSMatName[m]));
141: } else PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[m]));
143: for (i=0;i<rows;i++) {
144: 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: 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: PetscCall(PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v));
156: #endif
157: v += ds->ld;
158: #if defined(PETSC_USE_COMPLEX)
159: if (m==DS_MAT_T) vr += ds->ld;
160: #endif
161: }
162: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
163: }
164: PetscCall(MatDenseRestoreArrayRead(ds->omat[m],&M));
166: if (format == PETSC_VIEWER_ASCII_MATLAB) PetscCall(PetscViewerASCIIPrintf(viewer,"];\n"));
167: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
168: PetscCall(PetscViewerFlush(viewer));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: PetscErrorCode DSSortEigenvalues_Private(DS ds,PetscScalar *wr,PetscScalar *wi,PetscInt *perm,PetscBool isghiep)
173: {
174: PetscScalar re,im,wi0;
175: PetscInt n,i,j,result,tmp1,tmp2=0,d=1;
177: PetscFunctionBegin;
178: n = ds->t; /* sort only first t pairs if truncated */
179: /* insertion sort */
180: i=ds->l+1;
181: #if !defined(PETSC_USE_COMPLEX)
182: 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: for (;i<n;i+=d) {
187: re = wr[perm[i]];
188: if (wi) im = wi[perm[i]];
189: else im = 0.0;
190: tmp1 = perm[i];
191: #if !defined(PETSC_USE_COMPLEX)
192: 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: j = i-1;
199: if (wi) wi0 = wi[perm[j]];
200: else wi0 = 0.0;
201: PetscCall(SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result));
202: while (result<0 && j>=ds->l) {
203: perm[j+d] = perm[j];
204: j--;
205: #if !defined(PETSC_USE_COMPLEX)
206: if (wi && wi[perm[j+1]]!=0)
207: #else
208: if (isghiep && PetscImaginaryPart(wr[perm[j+1]])!=0)
209: #endif
210: { perm[j+d] = perm[j]; j--; }
212: if (j>=ds->l) {
213: if (wi) wi0 = wi[perm[j]];
214: else wi0 = 0.0;
215: PetscCall(SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result));
216: }
217: }
218: perm[j+1] = tmp1;
219: if (d==2) perm[j+2] = tmp2;
220: }
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: PetscErrorCode DSSortEigenvaluesReal_Private(DS ds,PetscReal *eig,PetscInt *perm)
225: {
226: PetscScalar re;
227: PetscInt i,j,result,tmp,l,n;
229: PetscFunctionBegin;
230: n = ds->t; /* sort only first t pairs if truncated */
231: l = ds->l;
232: /* insertion sort */
233: for (i=l+1;i<n;i++) {
234: re = eig[perm[i]];
235: j = i-1;
236: PetscCall(SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result));
237: while (result<0 && j>=l) {
238: tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
239: if (j>=l) PetscCall(SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result));
240: }
241: }
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*
246: Permute comumns [istart..iend-1] of [mat] according to perm. Columns have length n
247: */
248: PetscErrorCode DSPermuteColumns_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat,PetscInt *perm)
249: {
250: PetscInt i,j,k,p,ld;
251: PetscScalar *Q,rtmp;
253: PetscFunctionBegin;
254: ld = ds->ld;
255: PetscCall(MatDenseGetArray(ds->omat[mat],&Q));
256: for (i=istart;i<iend;i++) {
257: p = perm[i];
258: if (p != i) {
259: j = i + 1;
260: while (perm[j] != i) j++;
261: perm[j] = p; perm[i] = i;
262: /* swap columns i and j */
263: for (k=0;k<n;k++) {
264: rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
265: }
266: }
267: }
268: PetscCall(MatDenseRestoreArray(ds->omat[mat],&Q));
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: /*
273: The same as DSPermuteColumns_Private but for two matrices [mat1] and [mat2]
274: */
275: PetscErrorCode DSPermuteColumnsTwo_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat1,DSMatType mat2,PetscInt *perm)
276: {
277: PetscInt i,j,k,p,ld;
278: PetscScalar *Q,*Z,rtmp,rtmp2;
280: PetscFunctionBegin;
281: ld = ds->ld;
282: PetscCall(MatDenseGetArray(ds->omat[mat1],&Q));
283: PetscCall(MatDenseGetArray(ds->omat[mat2],&Z));
284: for (i=istart;i<iend;i++) {
285: p = perm[i];
286: if (p != i) {
287: j = i + 1;
288: while (perm[j] != i) j++;
289: perm[j] = p; perm[i] = i;
290: /* swap columns i and j */
291: for (k=0;k<n;k++) {
292: rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
293: rtmp2 = Z[k+p*ld]; Z[k+p*ld] = Z[k+i*ld]; Z[k+i*ld] = rtmp2;
294: }
295: }
296: }
297: PetscCall(MatDenseRestoreArray(ds->omat[mat1],&Q));
298: PetscCall(MatDenseRestoreArray(ds->omat[mat2],&Z));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: /*
303: Permute rows [istart..iend-1] of [mat] according to perm. Rows have length m
304: */
305: PetscErrorCode DSPermuteRows_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt m,DSMatType mat,PetscInt *perm)
306: {
307: PetscInt i,j,k,p,ld;
308: PetscScalar *Q,rtmp;
310: PetscFunctionBegin;
311: ld = ds->ld;
312: PetscCall(MatDenseGetArray(ds->omat[mat],&Q));
313: for (i=istart;i<iend;i++) {
314: p = perm[i];
315: if (p != i) {
316: j = i + 1;
317: while (perm[j] != i) j++;
318: perm[j] = p; perm[i] = i;
319: /* swap rows i and j */
320: for (k=0;k<m;k++) {
321: rtmp = Q[p+k*ld]; Q[p+k*ld] = Q[i+k*ld]; Q[i+k*ld] = rtmp;
322: }
323: }
324: }
325: PetscCall(MatDenseRestoreArray(ds->omat[mat],&Q));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
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: PetscErrorCode DSPermuteBoth_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,PetscInt m,DSMatType mat1,DSMatType mat2,PetscInt *perm)
334: {
335: PetscInt i,j,k,p,ld;
336: PetscScalar *U,*V,rtmp;
338: PetscFunctionBegin;
339: ld = ds->ld;
340: PetscCall(MatDenseGetArray(ds->omat[mat1],&U));
341: PetscCall(MatDenseGetArray(ds->omat[mat2],&V));
342: for (i=istart;i<iend;i++) {
343: p = perm[i];
344: if (p != i) {
345: j = i + 1;
346: while (perm[j] != i) j++;
347: perm[j] = p; perm[i] = i;
348: /* swap columns i and j of U */
349: for (k=0;k<n;k++) {
350: 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: for (k=0;k<m;k++) {
354: rtmp = V[k+p*ld]; V[k+p*ld] = V[k+i*ld]; V[k+i*ld] = rtmp;
355: }
356: }
357: }
358: PetscCall(MatDenseRestoreArray(ds->omat[mat1],&U));
359: PetscCall(MatDenseRestoreArray(ds->omat[mat2],&V));
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: /*@
364: DSSetIdentity - Copy the identity (a diagonal matrix with ones) on the
365: active part of a matrix.
367: Logically Collective
369: Input Parameters:
370: + ds - the direct solver context
371: - mat - the matrix to modify
373: Level: intermediate
375: .seealso: DSGetMat()
376: @*/
377: PetscErrorCode DSSetIdentity(DS ds,DSMatType mat)
378: {
379: PetscScalar *A;
380: PetscInt i,ld,n,l;
382: PetscFunctionBegin;
385: DSCheckValidMat(ds,mat,2);
387: PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
388: PetscCall(DSGetLeadingDimension(ds,&ld));
389: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
390: PetscCall(MatDenseGetArray(ds->omat[mat],&A));
391: PetscCall(PetscArrayzero(A+l*ld,ld*(n-l)));
392: for (i=l;i<n;i++) A[i+i*ld] = 1.0;
393: PetscCall(MatDenseRestoreArray(ds->omat[mat],&A));
394: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: /*@
399: DSOrthogonalize - Orthogonalize the columns of a matrix.
401: Logically Collective
403: Input Parameters:
404: + ds - the direct solver context
405: . mat - a matrix
406: - cols - number of columns to orthogonalize (starting from column zero)
408: Output Parameter:
409: . lindcols - (optional) number of linearly independent columns of the matrix
411: Level: developer
413: .seealso: DSPseudoOrthogonalize()
414: @*/
415: PetscErrorCode DSOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscInt *lindcols)
416: {
417: PetscInt n,l,ld;
418: PetscBLASInt ld_,rA,cA,info,ltau,lw;
419: PetscScalar *A,*tau,*w,saux,dummy;
421: PetscFunctionBegin;
423: DSCheckAlloc(ds,1);
425: DSCheckValidMat(ds,mat,2);
428: PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
429: PetscCall(DSGetLeadingDimension(ds,&ld));
430: n = n - l;
431: PetscCheck(cols<=n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
432: if (n == 0 || cols == 0) PetscFunctionReturn(PETSC_SUCCESS);
434: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
435: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
436: PetscCall(MatDenseGetArray(ds->omat[mat],&A));
437: PetscCall(PetscBLASIntCast(PetscMin(cols,n),<au));
438: PetscCall(PetscBLASIntCast(ld,&ld_));
439: PetscCall(PetscBLASIntCast(n,&rA));
440: PetscCall(PetscBLASIntCast(cols,&cA));
441: lw = -1;
442: PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,A,&ld_,&dummy,&saux,&lw,&info));
443: SlepcCheckLapackInfo("geqrf",info);
444: lw = (PetscBLASInt)PetscRealPart(saux);
445: PetscCall(DSAllocateWork_Private(ds,lw+ltau,0,0));
446: tau = ds->work;
447: w = &tau[ltau];
448: PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,&A[ld*l+l],&ld_,tau,w,&lw,&info));
449: SlepcCheckLapackInfo("geqrf",info);
450: PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&rA,<au,<au,&A[ld*l+l],&ld_,tau,w,&lw,&info));
451: SlepcCheckLapackInfo("orgqr",info);
452: if (lindcols) *lindcols = ltau;
453: PetscCall(PetscFPTrapPop());
454: PetscCall(MatDenseRestoreArray(ds->omat[mat],&A));
455: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
456: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
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: 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: PetscInt tmp;
473: PetscBLASInt m, n, k, ldA, ldB, ldC;
474: const char *N = "N", *T = "C", *qA = N, *qB = N;
476: PetscFunctionBegin;
477: if ((rA == 0) || (cB == 0)) PetscFunctionReturn(PETSC_SUCCESS);
478: PetscAssertPointer(C,1);
479: PetscAssertPointer(A,5);
480: PetscAssertPointer(B,10);
481: PetscCall(PetscBLASIntCast(_ldA,&ldA));
482: PetscCall(PetscBLASIntCast(_ldB,&ldB));
483: PetscCall(PetscBLASIntCast(_ldC,&ldC));
485: /* Transpose if needed */
486: if (At) tmp = rA, rA = cA, cA = tmp, qA = T;
487: if (Bt) tmp = rB, rB = cB, cB = tmp, qB = T;
489: /* Check size */
490: PetscCheck(cA==rB,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Matrix dimensions do not match");
492: /* Do stub */
493: if ((rA == 1) && (cA == 1) && (cB == 1)) {
494: if (!At && !Bt) *C = *A * *B;
495: else if (At && !Bt) *C = PetscConj(*A) * *B;
496: else if (!At && Bt) *C = *A * PetscConj(*B);
497: else *C = PetscConj(*A) * PetscConj(*B);
498: m = n = k = 1;
499: } else {
500: PetscCall(PetscBLASIntCast(rA,&m));
501: PetscCall(PetscBLASIntCast(cB,&n));
502: PetscCall(PetscBLASIntCast(cA,&k));
503: PetscCallBLAS("BLASgemm",BLASgemm_(qA,qB,&m,&n,&k,&a,(PetscScalar*)A,&ldA,(PetscScalar*)B,&ldB,&b,C,&ldC));
504: }
506: PetscCall(PetscLogFlops(2.0*m*n*k));
507: PetscFunctionReturn(PETSC_SUCCESS);
508: }
510: /*@
511: DSPseudoOrthogonalize - Orthogonalize the columns of a matrix with Modified
512: Gram-Schmidt in an indefinite inner product space defined by a signature.
514: Logically Collective
516: Input Parameters:
517: + ds - the direct solver context
518: . mat - the matrix
519: . cols - number of columns to orthogonalize (starting from column zero)
520: - s - the signature that defines the inner product
522: Output Parameters:
523: + lindcols - (optional) linearly independent columns of the matrix
524: - ns - (optional) the new signature of the vectors
526: Note:
527: After the call the matrix satisfies A'*s*A = ns.
529: Level: developer
531: .seealso: DSOrthogonalize()
532: @*/
533: PetscErrorCode DSPseudoOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscReal s[],PetscInt *lindcols,PetscReal ns[])
534: {
535: PetscInt i,j,k,l,n,ld;
536: PetscBLASInt info,one=1,zero=0,rA_,ld_;
537: PetscScalar *A,*A_,*m,*h,nr0;
538: PetscReal nr_o,nr,nr_abs,*ns_,done=1.0;
540: PetscFunctionBegin;
542: DSCheckAlloc(ds,1);
544: DSCheckValidMat(ds,mat,2);
546: PetscAssertPointer(s,4);
547: PetscCall(DSGetDimensions(ds,&n,&l,NULL,NULL));
548: PetscCall(DSGetLeadingDimension(ds,&ld));
549: n = n - l;
550: PetscCheck(cols<=n,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
551: if (n == 0 || cols == 0) PetscFunctionReturn(PETSC_SUCCESS);
552: PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
553: PetscCall(PetscBLASIntCast(n,&rA_));
554: PetscCall(PetscBLASIntCast(ld,&ld_));
555: PetscCall(MatDenseGetArray(ds->omat[mat],&A_));
556: A = &A_[ld*l+l];
557: PetscCall(DSAllocateWork_Private(ds,n+cols,ns?0:cols,0));
558: m = ds->work;
559: h = &m[n];
560: ns_ = ns ? ns : ds->rwork;
561: for (i=0; i<cols; i++) {
562: /* m <- diag(s)*A[i] */
563: for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
564: /* nr_o <- mynorm(A[i]'*m), mynorm(x) = sign(x)*sqrt(|x|) */
565: PetscCall(SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
566: nr = nr_o = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
567: for (j=0; j<3 && i>0; j++) {
568: /* h <- A[0:i-1]'*m */
569: PetscCall(SlepcMatDenseMult(h,i,0.0,1.0,A,ld,n,i,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
570: /* h <- diag(ns)*h */
571: for (k=0; k<i; k++) h[k] *= ns_[k];
572: /* A[i] <- A[i] - A[0:i-1]*h */
573: PetscCall(SlepcMatDenseMult(&A[ld*i],ld,1.0,-1.0,A,ld,n,i,PETSC_FALSE,h,i,i,1,PETSC_FALSE));
574: /* m <- diag(s)*A[i] */
575: for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
576: /* nr_o <- mynorm(A[i]'*m) */
577: PetscCall(SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE));
578: nr = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
579: PetscCheck(PetscAbs(nr)>PETSC_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Linear dependency detected");
580: if (PetscAbs(nr) > 0.7*PetscAbs(nr_o)) break;
581: nr_o = nr;
582: }
583: ns_[i] = PetscSign(nr);
584: /* A[i] <- A[i]/|nr| */
585: nr_abs = PetscAbs(nr);
586: PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&nr_abs,&done,&rA_,&one,A+i*ld,&ld_,&info));
587: SlepcCheckLapackInfo("lascl",info);
588: }
589: PetscCall(MatDenseRestoreArray(ds->omat[mat],&A_));
590: PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
591: PetscCall(PetscObjectStateIncrease((PetscObject)ds));
592: if (lindcols) *lindcols = cols;
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }