Actual source code: bvlapack.c

slepc-3.14.1 2020-12-08
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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:    BV private kernels that use the LAPACK
 12: */

 14: #include <slepc/private/bvimpl.h>
 15: #include <slepcblaslapack.h>

 17: /*
 18:     Reduction operation to compute sqrt(x**2+y**2) when normalizing vectors
 19: */
 20: SLEPC_EXTERN void MPIAPI SlepcPythag(void *in,void *inout,PetscMPIInt *len,MPI_Datatype *datatype)
 21: {
 22:   PetscBLASInt i,n=*len;
 23:   PetscReal    *x = (PetscReal*)in,*y = (PetscReal*)inout;

 26:   if (PetscUnlikely(*datatype!=MPIU_REAL)) {
 27:     (*PetscErrorPrintf)("Only implemented for MPIU_REAL data type");
 28:     MPI_Abort(MPI_COMM_WORLD,1);
 29:   }
 30:   for (i=0;i<n;i++) y[i] = SlepcAbs(x[i],y[i]);
 31:   PetscFunctionReturnVoid();
 32: }

 34: /*
 35:     Compute ||A|| for an mxn matrix
 36: */
 37: PetscErrorCode BVNorm_LAPACK_Private(BV bv,PetscInt m_,PetscInt n_,const PetscScalar *A,NormType type,PetscReal *nrm,PetscBool mpi)
 38: {
 40:   PetscBLASInt   m,n,i,j;
 41:   PetscMPIInt    len;
 42:   PetscReal      lnrm,*rwork=NULL,*rwork2=NULL;

 45:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 46:   PetscBLASIntCast(m_,&m);
 47:   PetscBLASIntCast(n_,&n);
 48:   if (type==NORM_FROBENIUS || type==NORM_2) {
 49:     lnrm = LAPACKlange_("F",&m,&n,(PetscScalar*)A,&m,rwork);
 50:     if (mpi) {
 51:       MPIU_Allreduce(&lnrm,nrm,1,MPIU_REAL,MPIU_LAPY2,PetscObjectComm((PetscObject)bv));
 52:     } else *nrm = lnrm;
 53:     PetscLogFlops(2.0*m*n);
 54:   } else if (type==NORM_1) {
 55:     if (mpi) {
 56:       BVAllocateWork_Private(bv,2*n_);
 57:       rwork = (PetscReal*)bv->work;
 58:       rwork2 = rwork+n_;
 59:       PetscArrayzero(rwork,n_);
 60:       PetscArrayzero(rwork2,n_);
 61:       for (j=0;j<n_;j++) {
 62:         for (i=0;i<m_;i++) {
 63:           rwork[j] += PetscAbsScalar(A[i+j*m_]);
 64:         }
 65:       }
 66:       PetscMPIIntCast(n_,&len);
 67:       MPIU_Allreduce(rwork,rwork2,len,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)bv));
 68:       *nrm = 0.0;
 69:       for (j=0;j<n_;j++) if (rwork2[j] > *nrm) *nrm = rwork2[j];
 70:     } else {
 71:       *nrm = LAPACKlange_("O",&m,&n,(PetscScalar*)A,&m,rwork);
 72:     }
 73:     PetscLogFlops(1.0*m*n);
 74:   } else if (type==NORM_INFINITY) {
 75:     BVAllocateWork_Private(bv,m_);
 76:     rwork = (PetscReal*)bv->work;
 77:     lnrm = LAPACKlange_("I",&m,&n,(PetscScalar*)A,&m,rwork);
 78:     if (mpi) {
 79:       MPIU_Allreduce(&lnrm,nrm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)bv));
 80:     } else *nrm = lnrm;
 81:     PetscLogFlops(1.0*m*n);
 82:   }
 83:   PetscFPTrapPop();
 84:   return(0);
 85: }

 87: /*
 88:     Normalize the columns of an mxn matrix A
 89: */
 90: PetscErrorCode BVNormalize_LAPACK_Private(BV bv,PetscInt m_,PetscInt n_,const PetscScalar *A,PetscScalar *eigi,PetscBool mpi)
 91: {
 93:   PetscBLASInt   m,j,k,info,zero=0;
 94:   PetscMPIInt    len;
 95:   PetscReal      *norms,*rwork=NULL,*rwork2=NULL,done=1.0;

 98:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 99:   PetscBLASIntCast(m_,&m);
100:   BVAllocateWork_Private(bv,2*n_);
101:   rwork = (PetscReal*)bv->work;
102:   rwork2 = rwork+n_;
103:   /* compute local norms */
104:   for (j=0;j<n_;j++) {
105:     k = 1;
106: #if !defined(PETSC_USE_COMPLEX)
107:     if (eigi && eigi[j] != 0.0) k = 2;
108: #endif
109:     rwork[j] = LAPACKlange_("F",&m,&k,(PetscScalar*)(A+j*m_),&m,rwork2);
110:     if (k==2) { rwork[j+1] = rwork[j]; j++; }
111:   }
112:   /* reduction to get global norms */
113:   if (mpi) {
114:     PetscMPIIntCast(n_,&len);
115:     PetscArrayzero(rwork2,n_);
116:     MPIU_Allreduce(rwork,rwork2,len,MPIU_REAL,MPIU_LAPY2,PetscObjectComm((PetscObject)bv));
117:     norms = rwork2;
118:   } else norms = rwork;
119:   /* scale columns */
120:   for (j=0;j<n_;j++) {
121:     k = 1;
122: #if !defined(PETSC_USE_COMPLEX)
123:     if (eigi && eigi[j] != 0.0) k = 2;
124: #endif
125:     PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,norms+j,&done,&m,&k,(PetscScalar*)(A+j*m_),&m,&info));
126:     SlepcCheckLapackInfo("lascl",info);
127:     if (k==2) j++;
128:   }
129:   PetscLogFlops(3.0*m*n_);
130:   PetscFPTrapPop();
131:   return(0);
132: }

134: /*
135:    Compute the upper Cholesky factor in R and its inverse in S.
136:    If S == R then the inverse overwrites the Cholesky factor.
137:  */
138: PetscErrorCode BVMatCholInv_LAPACK_Private(BV bv,Mat R,Mat S)
139: {
141:   PetscInt       i,k,l,n,m,ld,lds;
142:   PetscScalar    *pR,*pS;
143:   PetscBLASInt   info,n_ = 0,m_ = 0,ld_,lds_;

146:   l = bv->l;
147:   k = bv->k;
148:   MatGetSize(R,&m,NULL);
149:   n = k-l;
150:   PetscBLASIntCast(m,&m_);
151:   PetscBLASIntCast(n,&n_);
152:   ld  = m;
153:   ld_ = m_;
154:   MatDenseGetArray(R,&pR);

156:   if (S==R) {
157:     BVAllocateWork_Private(bv,m*k);
158:     pS = bv->work;
159:     lds = ld;
160:     lds_ = ld_;
161:   } else {
162:     MatDenseGetArray(S,&pS);
163:     MatGetSize(S,&lds,NULL);
164:     PetscBLASIntCast(lds,&lds_);
165:   }

167:   /* save a copy of matrix in S */
168:   for (i=l;i<k;i++) {
169:     PetscArraycpy(pS+i*lds+l,pR+i*ld+l,n);
170:   }

172:   /* compute upper Cholesky factor in R */
173:   PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&n_,pR+l*ld+l,&ld_,&info));
174:   PetscLogFlops((1.0*n*n*n)/3.0);

176:   if (info) {  /* LAPACKpotrf failed, retry on diagonally perturbed matrix */
177:     for (i=l;i<k;i++) {
178:       PetscArraycpy(pR+i*ld+l,pS+i*lds+l,n);
179:       pR[i+i*ld] += 50.0*PETSC_MACHINE_EPSILON;
180:     }
181:     PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("U",&n_,pR+l*ld+l,&ld_,&info));
182:     SlepcCheckLapackInfo("potrf",info);
183:     PetscLogFlops((1.0*n*n*n)/3.0);
184:   }

186:   /* compute S = inv(R) */
187:   if (S==R) {
188:     PetscStackCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&n_,pR+l*ld+l,&ld_,&info));
189:   } else {
190:     PetscArrayzero(pS+l*lds,(k-l)*k);
191:     for (i=l;i<k;i++) {
192:       PetscArraycpy(pS+i*lds+l,pR+i*ld+l,n);
193:     }
194:     PetscStackCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&n_,pS+l*lds+l,&lds_,&info));
195:   }
196:   SlepcCheckLapackInfo("trtri",info);
197:   PetscLogFlops(0.33*n*n*n);

199:   /* Zero out entries below the diagonal */
200:   for (i=l;i<k-1;i++) {
201:     PetscArrayzero(pR+i*ld+i+1,(k-i-1));
202:     if (S!=R) { PetscArrayzero(pS+i*lds+i+1,(k-i-1)); }
203:   }
204:   MatDenseRestoreArray(R,&pR);
205:   if (S!=R) { MatDenseRestoreArray(S,&pS); }
206:   return(0);
207: }

209: /*
210:    Compute the inverse of an upper triangular matrix R, store it in S.
211:    If S == R then the inverse overwrites R.
212:  */
213: PetscErrorCode BVMatTriInv_LAPACK_Private(BV bv,Mat R,Mat S)
214: {
216:   PetscInt       i,k,l,n,m,ld,lds;
217:   PetscScalar    *pR,*pS;
218:   PetscBLASInt   info,n_,m_ = 0,ld_,lds_;

221:   l = bv->l;
222:   k = bv->k;
223:   MatGetSize(R,&m,NULL);
224:   n = k-l;
225:   PetscBLASIntCast(m,&m_);
226:   PetscBLASIntCast(n,&n_);
227:   ld  = m;
228:   ld_ = m_;
229:   MatDenseGetArray(R,&pR);

231:   if (S==R) {
232:     BVAllocateWork_Private(bv,m*k);
233:     pS = bv->work;
234:     lds = ld;
235:     lds_ = ld_;
236:   } else {
237:     MatDenseGetArray(S,&pS);
238:     MatGetSize(S,&lds,NULL);
239:     PetscBLASIntCast(lds,&lds_);
240:   }

242:   /* compute S = inv(R) */
243:   if (S==R) {
244:     PetscStackCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&n_,pR+l*ld+l,&ld_,&info));
245:   } else {
246:     PetscArrayzero(pS+l*lds,(k-l)*k);
247:     for (i=l;i<k;i++) {
248:       PetscArraycpy(pS+i*lds+l,pR+i*ld+l,n);
249:     }
250:     PetscStackCallBLAS("LAPACKtrtri",LAPACKtrtri_("U","N",&n_,pS+l*lds+l,&lds_,&info));
251:   }
252:   SlepcCheckLapackInfo("trtri",info);
253:   PetscLogFlops(0.33*n*n*n);

255:   MatDenseRestoreArray(R,&pR);
256:   if (S!=R) { MatDenseRestoreArray(S,&pS); }
257:   return(0);
258: }

260: /*
261:    Compute the matrix to be used for post-multiplying the basis in the SVQB
262:    block orthogonalization method.
263:    On input R = V'*V, on output S = D*U*Lambda^{-1/2} where (U,Lambda) is
264:    the eigendecomposition of D*R*D with D=diag(R)^{-1/2}.
265:    If S == R then the result overwrites R.
266:  */
267: PetscErrorCode BVMatSVQB_LAPACK_Private(BV bv,Mat R,Mat S)
268: {
270:   PetscInt       i,j,k,l,n,m,ld,lds;
271:   PetscScalar    *pR,*pS,*D,*work,a;
272:   PetscReal      *eig,dummy;
273:   PetscBLASInt   info,lwork,n_,m_ = 0,ld_,lds_;
274: #if defined(PETSC_USE_COMPLEX)
275:   PetscReal      *rwork,rdummy;
276: #endif

279:   l = bv->l;
280:   k = bv->k;
281:   MatGetSize(R,&m,NULL);
282:   n = k-l;
283:   PetscBLASIntCast(m,&m_);
284:   PetscBLASIntCast(n,&n_);
285:   ld  = m;
286:   ld_ = m_;
287:   MatDenseGetArray(R,&pR);

289:   if (S==R) {
290:     pS = pR;
291:     lds = ld;
292:     lds_ = ld_;
293:   } else {
294:     MatDenseGetArray(S,&pS);
295:     MatGetSize(S,&lds,NULL);
296:     PetscBLASIntCast(lds,&lds_);
297:   }

299:   /* workspace query and memory allocation */
300:   lwork = -1;
301: #if defined(PETSC_USE_COMPLEX)
302:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n_,pS,&lds_,&dummy,&a,&lwork,&rdummy,&info));
303:   PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
304:   PetscMalloc4(n,&eig,n,&D,lwork,&work,PetscMax(1,3*n-2),&rwork);
305: #else
306:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n_,pS,&lds_,&dummy,&a,&lwork,&info));
307:   PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork);
308:   PetscMalloc3(n,&eig,n,&D,lwork,&work);
309: #endif

311:   /* copy and scale matrix */
312:   for (i=l;i<k;i++) D[i-l] = 1.0/PetscSqrtReal(PetscRealPart(pR[i+i*ld]));
313:   for (i=l;i<k;i++) for (j=l;j<k;j++) pS[i+j*lds] = pR[i+j*ld]*D[i-l];
314:   for (j=l;j<k;j++) for (i=l;i<k;i++) pS[i+j*lds] *= D[j-l];

316:   /* compute eigendecomposition */
317: #if defined(PETSC_USE_COMPLEX)
318:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n_,pS+l*lds+l,&lds_,eig,work,&lwork,rwork,&info));
319: #else
320:   PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&n_,pS+l*lds+l,&lds_,eig,work,&lwork,&info));
321: #endif
322:   SlepcCheckLapackInfo("syev",info);

324:   if (S!=R) {   /* R = U' */
325:     for (i=l;i<k;i++) for (j=l;j<k;j++) pR[i+j*ld] = pS[j+i*lds];
326:   }

328:   /* compute S = D*U*Lambda^{-1/2} */
329:   for (i=l;i<k;i++) for (j=l;j<k;j++) pS[i+j*lds] *= D[i-l];
330:   for (j=l;j<k;j++) for (i=l;i<k;i++) pS[i+j*lds] /= PetscSqrtReal(eig[j-l]);

332:   if (S!=R) {   /* compute R = inv(S) = Lambda^{1/2}*U'/D */
333:     for (i=l;i<k;i++) for (j=l;j<k;j++) pR[i+j*ld] *= PetscSqrtReal(eig[i-l]);
334:     for (j=l;j<k;j++) for (i=l;i<k;i++) pR[i+j*ld] /= D[j-l];
335:   }

337: #if defined(PETSC_USE_COMPLEX)
338:   PetscFree4(eig,D,work,rwork);
339: #else
340:   PetscFree3(eig,D,work);
341: #endif
342:   PetscLogFlops(9.0*n*n*n);

344:   MatDenseRestoreArray(R,&pR);
345:   if (S!=R) { MatDenseRestoreArray(S,&pS); }
346:   return(0);
347: }

349: /*
350:     QR factorization of an mxn matrix via parallel TSQR
351: */
352: PetscErrorCode BVOrthogonalize_LAPACK_TSQR(BV bv,PetscInt m_,PetscInt n_,PetscScalar *Q,PetscScalar *R,PetscInt ldr)
353: {
355:   PetscInt       level,plevel,nlevels,powtwo,lda,worklen;
356:   PetscBLASInt   m,n,i,j,k,l,s = 0,nb,sz,lwork,info;
357:   PetscScalar    *tau,*work,*A=NULL,*QQ=NULL,*Qhalf,*C=NULL,one=1.0,zero=0.0;
358:   PetscMPIInt    rank,size,count,stride;
359:   MPI_Datatype   tmat;

362:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
363:   PetscBLASIntCast(m_,&m);
364:   PetscBLASIntCast(n_,&n);
365:   k  = PetscMin(m,n);
366:   nb = 16;
367:   lda = 2*n;
368:   MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size);
369:   MPI_Comm_rank(PetscObjectComm((PetscObject)bv),&rank);
370:   nlevels = (PetscInt)PetscCeilReal(PetscLog2Real((PetscReal)size));
371:   powtwo  = PetscPowInt(2,(PetscInt)PetscFloorReal(PetscLog2Real((PetscReal)size)));
372:   worklen = n+n*nb;
373:   if (nlevels) worklen += n*lda+n*lda*nlevels+n*lda;
374:   BVAllocateWork_Private(bv,worklen);
375:   tau  = bv->work;
376:   work = bv->work+n;
377:   PetscBLASIntCast(n*nb,&lwork);
378:   if (nlevels) {
379:     A  = bv->work+n+n*nb;
380:     QQ = bv->work+n+n*nb+n*lda;
381:     C  = bv->work+n+n*nb+n*lda+n*lda*nlevels;
382:   }

384:   /* Compute QR */
385:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,Q,&m,tau,work,&lwork,&info));
386:   SlepcCheckLapackInfo("geqrf",info);

388:   /* Extract R */
389:   if (R || nlevels) {
390:     for (j=0;j<n;j++) {
391:       for (i=0;i<=PetscMin(j,m-1);i++) {
392:         if (nlevels) A[i+j*lda] = Q[i+j*m];
393:         else R[i+j*ldr] = Q[i+j*m];
394:       }
395:       for (i=PetscMin(j,m-1)+1;i<n;i++) {
396:         if (nlevels) A[i+j*lda] = 0.0;
397:         else R[i+j*ldr] = 0.0;
398:       }
399:     }
400:   }

402:   /* Compute orthogonal matrix in Q */
403:   PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&m,&k,&k,Q,&m,tau,work,&lwork,&info));
404:   SlepcCheckLapackInfo("orgqr",info);

406:   if (nlevels) {

408:     PetscMPIIntCast(n,&count);
409:     PetscMPIIntCast(lda,&stride);
410:     PetscBLASIntCast(lda,&l);
411:     MPI_Type_vector(count,count,stride,MPIU_SCALAR,&tmat);
412:     MPI_Type_commit(&tmat);

414:     for (level=nlevels;level>=1;level--) {

416:       plevel = PetscPowInt(2,level);
417:       PetscBLASIntCast(plevel*PetscFloorReal(rank/(PetscReal)plevel)+(rank+PetscPowInt(2,level-1))%plevel,&s);

419:       /* Stack triangular matrices */
420:       if (rank<s && s<size) {  /* send top part, receive bottom part */
421:         MPI_Sendrecv(A,1,tmat,s,111,A+n,1,tmat,s,111,PetscObjectComm((PetscObject)bv),MPI_STATUS_IGNORE);
422:       } else if (s<size) {  /* copy top to bottom, receive top part */
423:         MPI_Sendrecv(A,1,tmat,rank,111,A+n,1,tmat,rank,111,PetscObjectComm((PetscObject)bv),MPI_STATUS_IGNORE);
424:         MPI_Sendrecv(A+n,1,tmat,s,111,A,1,tmat,s,111,PetscObjectComm((PetscObject)bv),MPI_STATUS_IGNORE);
425:       }
426:       if (level<nlevels && size!=powtwo) {  /* for cases when size is not a power of 2 */
427:         if (rank<size-powtwo) {  /* send bottom part */
428:           MPI_Send(A+n,1,tmat,rank+powtwo,111,PetscObjectComm((PetscObject)bv));
429:         } else if (rank>=powtwo) {  /* receive bottom part */
430:           MPI_Recv(A+n,1,tmat,rank-powtwo,111,PetscObjectComm((PetscObject)bv),MPI_STATUS_IGNORE);
431:         }
432:       }
433:       /* Compute QR and build orthogonal matrix */
434:       if (level<nlevels || (level==nlevels && s<size)) {
435:         PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&l,&n,A,&l,tau,work,&lwork,&info));
436:         SlepcCheckLapackInfo("geqrf",info);
437:         PetscArraycpy(QQ+(level-1)*n*lda,A,n*lda);
438:         PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&l,&n,&n,QQ+(level-1)*n*lda,&l,tau,work,&lwork,&info));
439:         SlepcCheckLapackInfo("orgqr",info);
440:         for (j=0;j<n;j++) {
441:           for (i=j+1;i<n;i++) A[i+j*lda] = 0.0;
442:         }
443:       } else if (level==nlevels) {  /* only one triangular matrix, set Q=I */
444:         PetscArrayzero(QQ+(level-1)*n*lda,n*lda);
445:         for (j=0;j<n;j++) QQ[j+j*lda+(level-1)*n*lda] = 1.0;
446:       }
447:     }

449:     /* Extract R */
450:     if (R) {
451:       for (j=0;j<n;j++) {
452:         for (i=0;i<=j;i++) R[i+j*ldr] = A[i+j*lda];
453:         for (i=j+1;i<n;i++) R[i+j*ldr] = 0.0;
454:       }
455:     }

457:     /* Accumulate orthogonal matrices */
458:     for (level=1;level<=nlevels;level++) {
459:       plevel = PetscPowInt(2,level);
460:       PetscBLASIntCast(plevel*PetscFloorReal(rank/(PetscReal)plevel)+(rank+PetscPowInt(2,level-1))%plevel,&s);
461:       Qhalf = (rank<s)? QQ+(level-1)*n*lda: QQ+(level-1)*n*lda+n;
462:       if (level<nlevels) {
463:         PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&n,&one,QQ+level*n*lda,&l,Qhalf,&l,&zero,C,&l));
464:         PetscArraycpy(QQ+level*n*lda,C,n*lda);
465:       } else {
466:         for (i=0;i<m/l;i++) {
467:           PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&l,&n,&n,&one,Q+i*l,&m,Qhalf,&l,&zero,C,&l));
468:           for (j=0;j<n;j++) { PetscArraycpy(Q+i*l+j*m,C+j*l,l); }
469:         }
470:         sz = m%l;
471:         if (sz) {
472:           PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&sz,&n,&n,&one,Q+(m/l)*l,&m,Qhalf,&l,&zero,C,&l));
473:           for (j=0;j<n;j++) { PetscArraycpy(Q+(m/l)*l+j*m,C+j*l,sz); }
474:         }
475:       }
476:     }

478:     MPI_Type_free(&tmat);
479:   }

481:   PetscLogFlops(3.0*m*n*n);
482:   PetscFPTrapPop();
483:   return(0);
484: }

486: /*
487:     Reduction operation to compute [~,Rout]=qr([Rin1;Rin2]) in the TSQR algorithm;
488:     all matrices are upper triangular stored in packed format
489: */
490: SLEPC_EXTERN void MPIAPI SlepcGivensPacked(void *in,void *inout,PetscMPIInt *len,MPI_Datatype *datatype)
491: {
493:   PetscBLASInt   n,i,j,k,one=1;
494:   PetscMPIInt    tsize;
495:   PetscScalar    v,s,*R2=(PetscScalar*)in,*R1=(PetscScalar*)inout;
496:   PetscReal      c;

499:   MPI_Type_size(*datatype,&tsize);CHKERRABORT(PETSC_COMM_SELF,ierr);  /* we assume len=1 */
500:   tsize /= sizeof(PetscScalar);
501:   n = (-1+(PetscBLASInt)PetscSqrtReal(1+8*tsize))/2;
502:   for (j=0;j<n;j++) {
503:     for (i=0;i<=j;i++) {
504:       LAPACKlartg_(R1+(2*n-j-1)*j/2+j,R2+(2*n-i-1)*i/2+j,&c,&s,&v);
505:       R1[(2*n-j-1)*j/2+j] = v;
506:       k = n-j-1;
507:       if (k) BLASrot_(&k,R1+(2*n-j-1)*j/2+j+1,&one,R2+(2*n-i-1)*i/2+j+1,&one,&c,&s);
508:     }
509:   }
510:   PetscFunctionReturnVoid();
511: }

513: /*
514:     Computes the R factor of the QR factorization of an mxn matrix via parallel TSQR
515: */
516: PetscErrorCode BVOrthogonalize_LAPACK_TSQR_OnlyR(BV bv,PetscInt m_,PetscInt n_,PetscScalar *Q,PetscScalar *R,PetscInt ldr)
517: {
519:   PetscInt       worklen;
520:   PetscBLASInt   m,n,i,j,s,nb,lwork,info;
521:   PetscScalar    *tau,*work,*A=NULL,*R1=NULL,*R2=NULL;
522:   PetscMPIInt    size,count;
523:   MPI_Datatype   tmat;

526:   PetscFPTrapPush(PETSC_FP_TRAP_OFF);
527:   PetscBLASIntCast(m_,&m);
528:   PetscBLASIntCast(n_,&n);
529:   nb = 16;
530:   s  = n+n*(n-1)/2;  /* length of packed triangular matrix */
531:   MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size);
532:   worklen = n+n*nb+2*s+m*n;
533:   BVAllocateWork_Private(bv,worklen);
534:   tau  = bv->work;
535:   work = bv->work+n;
536:   R1   = bv->work+n+n*nb;
537:   R2   = bv->work+n+n*nb+s;
538:   A    = bv->work+n+n*nb+2*s;
539:   PetscBLASIntCast(n*nb,&lwork);
540:   PetscArraycpy(A,Q,m*n);

542:   /* Compute QR */
543:   PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&m,&n,A,&m,tau,work,&lwork,&info));
544:   SlepcCheckLapackInfo("geqrf",info);

546:   if (size==1) {
547:     /* Extract R */
548:     for (j=0;j<n;j++) {
549:       for (i=0;i<=PetscMin(j,m-1);i++) R[i+j*ldr] = A[i+j*m];
550:       for (i=PetscMin(j,m-1)+1;i<n;i++) R[i+j*ldr] = 0.0;
551:     }
552:   } else {
553:     /* Use MPI reduction operation to obtain global R */
554:     PetscMPIIntCast(s,&count);
555:     MPI_Type_contiguous(count,MPIU_SCALAR,&tmat);
556:     MPI_Type_commit(&tmat);
557:     for (i=0;i<n;i++) {
558:       for (j=i;j<n;j++) R1[(2*n-i-1)*i/2+j] = (i<m)?A[i+j*m]:0.0;
559:     }
560:     MPIU_Allreduce(R1,R2,1,tmat,MPIU_TSQR,PetscObjectComm((PetscObject)bv));
561:     for (i=0;i<n;i++) {
562:       for (j=0;j<i;j++) R[i+j*ldr] = 0.0;
563:       for (j=i;j<n;j++) R[i+j*ldr] = R2[(2*n-i-1)*i/2+j];
564:     }
565:     MPI_Type_free(&tmat);
566:   }

568:   PetscLogFlops(3.0*m*n*n);
569:   PetscFPTrapPop();
570:   return(0);
571: }