Actual source code: dspriv.c

slepc-main 2024-11-22
Report Typos and Errors
  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),&ltau));
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,&ltau,&ltau,&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: }