1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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 orthogonalization routines
12: */
14: #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
16: /*
17: BV_NormVecOrColumn - Compute the 2-norm of the working vector, irrespective of
18: whether it is in a column or not
19: */
20: PETSC_STATIC_INLINE PetscErrorCode BV_NormVecOrColumn(BV bv,PetscInt j,Vec v,PetscReal *nrm) 21: {
25: if (v) { BVNormVec(bv,v,NORM_2,nrm); }
26: else { BVNormColumn(bv,j,NORM_2,nrm); }
27: return(0);
28: }
30: /*
31: BVDotColumnInc - Same as BVDotColumn() but also including column j, which
32: is multiplied by itself
33: */
34: PETSC_STATIC_INLINE PetscErrorCode BVDotColumnInc(BV X,PetscInt j,PetscScalar *q) 35: {
37: PetscInt ksave;
38: Vec y;
41: PetscLogEventBegin(BV_DotVec,X,0,0,0);
42: ksave = X->k;
43: X->k = j+1;
44: BVGetColumn(X,j,&y);
45: (*X->ops->dotvec)(X,y,q);
46: BVRestoreColumn(X,j,&y);
47: X->k = ksave;
48: PetscLogEventEnd(BV_DotVec,X,0,0,0);
49: return(0);
50: }
52: /*
53: BVOrthogonalizeMGS1 - Compute one step of Modified Gram-Schmidt
54: */
55: static PetscErrorCode BVOrthogonalizeMGS1(BV bv,PetscInt j,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onrm,PetscReal *nrm) 56: {
58: PetscInt i;
59: PetscScalar dot;
60: Vec vi,z,w=v;
61: const PetscScalar *omega;
64: if (!v) { BVGetColumn(bv,j,&w); }
65: if (onrm) { BVNormVec(bv,w,NORM_2,onrm); }
66: z = w;
67: if (bv->indef) {
68: VecGetArrayRead(bv->omega,&omega);
69: }
70: for (i=-bv->nc;i<j;i++) {
71: if (which && i>=0 && !which[i]) continue;
72: BVGetColumn(bv,i,&vi);
73: /* h_i = ( v, v_i ) */
74: if (bv->matrix) {
75: BV_IPMatMult(bv,w);
76: z = bv->Bx;
77: }
78: VecDot(z,vi,&dot);
79: /* v <- v - h_i v_i */
80: BV_SetValue(bv,i,0,c,dot);
81: if (bv->indef) dot /= PetscRealPart(omega[bv->nc+i]);
82: VecAXPY(w,-dot,vi);
83: BVRestoreColumn(bv,i,&vi);
84: }
85: if (nrm) { BVNormVec(bv,w,NORM_2,nrm); }
86: if (!v) { BVRestoreColumn(bv,j,&w); }
87: BV_AddCoefficients(bv,j,h,c);
88: if (bv->indef) {
89: VecRestoreArrayRead(bv->omega,&omega);
90: }
91: return(0);
92: }
94: /*
95: BVOrthogonalizeCGS1 - Compute |v'| (estimated), |v| and one step of CGS with
96: only one global synchronization
97: */
98: static PetscErrorCode BVOrthogonalizeCGS1(BV bv,PetscInt j,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onorm,PetscReal *norm) 99: {
101: PetscReal sum,beta;
104: /* h = W^* v ; alpha = (v, v) */
105: bv->k = j;
106: if (onorm || norm) {
107: if (!v) {
108: BVDotColumnInc(bv,j,c);
109: BV_SquareRoot(bv,j,c,&beta);
110: } else {
111: BVDotVec(bv,v,c);
112: BVNormVec(bv,v,NORM_2,&beta);
113: }
114: } else {
115: if (!v) { BVDotColumn(bv,j,c); }
116: else { BVDotVec(bv,v,c); }
117: }
119: /* q = v - V h */
120: if (bv->indef) { BV_ApplySignature(bv,j,c,PETSC_TRUE); }
121: if (!v) { BVMultColumn(bv,-1.0,1.0,j,c); }
122: else { BVMultVec(bv,-1.0,1.0,v,c); }
123: if (bv->indef) { BV_ApplySignature(bv,j,c,PETSC_FALSE); }
125: /* compute |v| */
126: if (onorm) *onorm = beta;
128: if (norm) {
129: if (bv->indef) {
130: BV_NormVecOrColumn(bv,j,v,norm);
131: } else {
132: /* estimate |v'| from |v| */
133: BV_SquareSum(bv,j,c,&sum);
134: *norm = beta*beta-sum;
135: if (*norm <= 0.0) {
136: BV_NormVecOrColumn(bv,j,v,norm);
137: } else *norm = PetscSqrtReal(*norm);
138: }
139: }
140: BV_AddCoefficients(bv,j,h,c);
141: return(0);
142: }
144: #define BVOrthogonalizeGS1(a,b,c,d,e,f,g,h) ((bv->ops->gramschmidt)?(*bv->ops->gramschmidt):(mgs?BVOrthogonalizeMGS1:BVOrthogonalizeCGS1))(a,b,c,d,e,f,g,h)146: /*
147: BVOrthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt
149: j - the index of the column to orthogonalize (cannot use both j and v)
150: v - the vector to orthogonalize (cannot use both j and v)
151: which - logical array indicating selected columns (only used in MGS)
152: norm - (optional) norm of the vector after being orthogonalized
153: lindep - (optional) flag indicating possible linear dependence
154: */
155: static PetscErrorCode BVOrthogonalizeGS(BV bv,PetscInt j,Vec v,PetscBool *which,PetscReal *norm,PetscBool *lindep)156: {
158: PetscScalar *h,*c,*omega;
159: PetscReal onrm,nrm;
160: PetscInt k,l;
161: PetscBool mgs,dolindep,signature;
164: if (v) {
165: k = bv->k;
166: h = bv->h;
167: c = bv->c;
168: } else {
169: k = j;
170: h = NULL;
171: c = NULL;
172: }
174: mgs = (bv->orthog_type==BV_ORTHOG_MGS)? PETSC_TRUE: PETSC_FALSE;
176: /* if indefinite inner product, skip the computation of lindep */
177: if (bv->indef && lindep) *lindep = PETSC_FALSE;
178: dolindep = (!bv->indef && lindep)? PETSC_TRUE: PETSC_FALSE;
180: /* if indefinite and we are orthogonalizing a column, the norm must always be computed */
181: signature = (bv->indef && !v)? PETSC_TRUE: PETSC_FALSE;
183: BV_CleanCoefficients(bv,k,h);
185: switch (bv->orthog_ref) {
187: case BV_ORTHOG_REFINE_IFNEEDED:
188: BVOrthogonalizeGS1(bv,k,v,which,h,c,&onrm,&nrm);
189: /* repeat if ||q|| < eta ||h|| */
190: l = 1;
191: while (l<3 && nrm && PetscAbsReal(nrm) < bv->orthog_eta*PetscAbsReal(onrm)) {
192: l++;
193: if (mgs||bv->indef) onrm = nrm;
194: BVOrthogonalizeGS1(bv,k,v,which,h,c,(mgs||bv->indef)?NULL:&onrm,&nrm);
195: }
196: /* linear dependence check: criterion not satisfied in the last iteration */
197: if (dolindep) *lindep = PetscNot(nrm && PetscAbsReal(nrm) >= bv->orthog_eta*PetscAbsReal(onrm));
198: break;
200: case BV_ORTHOG_REFINE_NEVER:
201: BVOrthogonalizeGS1(bv,k,v,which,h,c,NULL,NULL);
202: /* compute ||v|| */
203: if (norm || dolindep || signature) {
204: BV_NormVecOrColumn(bv,k,v,&nrm);
205: }
206: /* linear dependence check: just test for exactly zero norm */
207: if (dolindep) *lindep = PetscNot(nrm);
208: break;
210: case BV_ORTHOG_REFINE_ALWAYS:
211: BVOrthogonalizeGS1(bv,k,v,which,h,c,NULL,NULL);
212: BVOrthogonalizeGS1(bv,k,v,which,h,c,dolindep?&onrm:NULL,(norm||dolindep||signature)?&nrm:NULL);
213: /* linear dependence check: criterion not satisfied in the second iteration */
214: if (dolindep) *lindep = PetscNot(nrm && PetscAbsReal(nrm) >= bv->orthog_eta*PetscAbsReal(onrm));
215: break;
216: }
217: if (signature) {
218: VecGetArray(bv->omega,&omega);
219: omega[bv->nc+k] = (nrm<0.0)? -1.0: 1.0;
220: VecRestoreArray(bv->omega,&omega);
221: }
222: if (norm) {
223: *norm = nrm;
224: if (!v) { /* store norm value next to the orthogonalization coefficients */
225: if (dolindep && *lindep) { BV_SetValue(bv,k,k,h,0.0); }
226: else { BV_SetValue(bv,k,k,h,nrm); }
227: }
228: }
229: return(0);
230: }
232: /*@
233: BVOrthogonalizeVec - Orthogonalize a given vector with respect to all
234: active columns.
236: Collective on BV238: Input Parameters:
239: + bv - the basis vectors context
240: - v - the vector
242: Output Parameters:
243: + H - (optional) coefficients computed during orthogonalization
244: . norm - (optional) norm of the vector after being orthogonalized
245: - lindep - (optional) flag indicating that refinement did not improve the quality
246: of orthogonalization
248: Notes:
249: This function is equivalent to BVOrthogonalizeColumn() but orthogonalizes
250: a vector as an argument rather than taking one of the BV columns. The
251: vector is orthogonalized against all active columns (k) and the constraints.
252: If H is given, it must have enough space to store k-l coefficients, where l
253: is the number of leading columns.
255: In the case of an indefinite inner product, the lindep parameter is not
256: computed (set to false).
258: Level: advanced
260: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVSetActiveColumns(), BVGetNumConstraints()
261: @*/
262: PetscErrorCode BVOrthogonalizeVec(BV bv,Vec v,PetscScalar *H,PetscReal *norm,PetscBool *lindep)263: {
265: PetscInt ksave,lsave;
271: BVCheckSizes(bv,1);
275: PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0);
276: ksave = bv->k;
277: lsave = bv->l;
278: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
279: BV_AllocateCoeffs(bv);
280: BV_AllocateSignature(bv);
281: BVOrthogonalizeGS(bv,0,v,NULL,norm,lindep);
282: bv->k = ksave;
283: bv->l = lsave;
284: if (H) { BV_StoreCoefficients(bv,bv->k,bv->h,H); }
285: PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0);
286: return(0);
287: }
289: /*@
290: BVOrthogonalizeColumn - Orthogonalize one of the column vectors with respect to
291: the previous ones.
293: Collective on BV295: Input Parameters:
296: + bv - the basis vectors context
297: - j - index of column to be orthogonalized
299: Output Parameters:
300: + H - (optional) coefficients computed during orthogonalization
301: . norm - (optional) norm of the vector after being orthogonalized
302: - lindep - (optional) flag indicating that refinement did not improve the quality
303: of orthogonalization
305: Notes:
306: This function applies an orthogonal projector to project vector V[j] onto
307: the orthogonal complement of the span of the columns of V[0..j-1],
308: where V[.] are the vectors of BV. The columns V[0..j-1] are assumed to be
309: mutually orthonormal.
311: Leading columns V[0..l-1] also participate in the orthogonalization, as well
312: as the constraints. If H is given, it must have enough space to store
313: j-l+1 coefficients (the last coefficient will contain the value norm, unless
314: the norm argument is NULL).
316: If a non-standard inner product has been specified with BVSetMatrix(),
317: then the vector is B-orthogonalized, using the non-standard inner product
318: defined by matrix B. The output vector satisfies V[j]'*B*V[0..j-1] = 0.
320: This routine does not normalize the resulting vector, see BVOrthonormalizeColumn().
322: In the case of an indefinite inner product, the lindep parameter is not
323: computed (set to false).
325: Level: advanced
327: .seealso: BVSetOrthogonalization(), BVSetMatrix(), BVSetActiveColumns(), BVOrthogonalize(), BVOrthogonalizeVec(), BVGetNumConstraints(), BVOrthonormalizeColumn()
328: @*/
329: PetscErrorCode BVOrthogonalizeColumn(BV bv,PetscInt j,PetscScalar *H,PetscReal *norm,PetscBool *lindep)330: {
332: PetscInt ksave,lsave;
338: BVCheckSizes(bv,1);
339: if (j<0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
340: if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,bv->m);
342: PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0);
343: ksave = bv->k;
344: lsave = bv->l;
345: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
346: if (!bv->buffer) { BVGetBufferVec(bv,&bv->buffer); }
347: BV_AllocateSignature(bv);
348: BVOrthogonalizeGS(bv,j,NULL,NULL,norm,lindep);
349: bv->k = ksave;
350: bv->l = lsave;
351: if (H) { BV_StoreCoefficients(bv,j,NULL,H); }
352: PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0);
353: PetscObjectStateIncrease((PetscObject)bv);
354: return(0);
355: }
357: /*@
358: BVOrthonormalizeColumn - Orthonormalize one of the column vectors with respect to
359: the previous ones. This is equivalent to a call to BVOrthogonalizeColumn()
360: followed by a call to BVScaleColumn() with the reciprocal of the norm.
362: Collective on BV364: Input Parameters:
365: + bv - the basis vectors context
366: . j - index of column to be orthonormalized
367: - replace - whether it is allowed to set the vector randomly
369: Output Parameters:
370: + norm - (optional) norm of the vector after orthogonalization and before normalization
371: - lindep - (optional) flag indicating that linear dependence was determined during
372: orthogonalization
374: Notes:
375: This function first orthogonalizes vector V[j] with respect to V[0..j-1],
376: where V[.] are the vectors of BV. A byproduct of this computation is norm,
377: the norm of the vector after orthogonalization. Secondly, it scales the
378: vector with 1/norm, so that the resulting vector has unit norm.
380: If after orthogonalization the vector V[j] is exactly zero, it cannot be normalized
381: because norm=0. In that case, it could be left as zero or replaced by a random
382: vector that is then orthonormalized. The latter is achieved by setting the
383: argument replace to TRUE. The vector will be replaced by a random vector also
384: if lindep was set to TRUE, even if the norm is not exaclty zero.
386: If the vector has been replaced by a random vector, the output arguments norm and
387: lindep will be set according to the orthogonalization of this new vector.
389: Level: advanced
391: .seealso: BVOrthogonalizeColumn(), BVScaleColumn()
392: @*/
393: PetscErrorCode BVOrthonormalizeColumn(BV bv,PetscInt j,PetscBool replace,PetscReal *norm,PetscBool *lindep)394: {
396: PetscScalar alpha;
397: PetscReal nrm;
398: PetscInt ksave,lsave;
399: PetscBool lndep;
405: BVCheckSizes(bv,1);
406: if (j<0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
407: if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,bv->m);
409: /* orthogonalize */
410: PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0);
411: ksave = bv->k;
412: lsave = bv->l;
413: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
414: if (!bv->buffer) { BVGetBufferVec(bv,&bv->buffer); }
415: BV_AllocateSignature(bv);
416: BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep);
417: if (replace && (nrm==0.0 || lndep)) {
418: PetscInfo(bv,"Vector was linearly dependent, generating a new random vector\n");
419: BVSetRandomColumn(bv,j);
420: BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep);
421: if (nrm==0.0 || lndep) { /* yet another attempt */
422: BVSetRandomColumn(bv,j);
423: BVOrthogonalizeGS(bv,j,NULL,NULL,&nrm,&lndep);
424: }
425: }
426: bv->k = ksave;
427: bv->l = lsave;
428: PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0);
430: /* scale */
431: if (nrm!=1.0 && nrm!=0.0) {
432: alpha = 1.0/nrm;
433: PetscLogEventBegin(BV_Scale,bv,0,0,0);
434: if (bv->n) {
435: (*bv->ops->scale)(bv,j,alpha);
436: }
437: PetscLogEventEnd(BV_Scale,bv,0,0,0);
438: }
439: if (norm) *norm = nrm;
440: if (lindep) *lindep = lndep;
441: PetscObjectStateIncrease((PetscObject)bv);
442: return(0);
443: }
445: /*@
446: BVOrthogonalizeSomeColumn - Orthogonalize one of the column vectors with
447: respect to some of the previous ones.
449: Collective on BV451: Input Parameters:
452: + bv - the basis vectors context
453: . j - index of column to be orthogonalized
454: - which - logical array indicating selected columns
456: Output Parameters:
457: + H - (optional) coefficients computed during orthogonalization
458: . norm - (optional) norm of the vector after being orthogonalized
459: - lindep - (optional) flag indicating that refinement did not improve the quality
460: of orthogonalization
462: Notes:
463: This function is similar to BVOrthogonalizeColumn(), but V[j] is
464: orthogonalized only against columns V[i] having which[i]=PETSC_TRUE.
465: The length of array which must be j at least.
467: The use of this operation is restricted to MGS orthogonalization type.
469: In the case of an indefinite inner product, the lindep parameter is not
470: computed (set to false).
472: Level: advanced
474: .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization()
475: @*/
476: PetscErrorCode BVOrthogonalizeSomeColumn(BV bv,PetscInt j,PetscBool *which,PetscScalar *H,PetscReal *norm,PetscBool *lindep)477: {
479: PetscInt ksave,lsave;
486: BVCheckSizes(bv,1);
487: if (j<0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
488: if (j>=bv->m) SETERRQ2(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but BV only has %D columns",j,bv->m);
489: if (bv->orthog_type!=BV_ORTHOG_MGS) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation only available for MGS orthogonalization");
491: PetscLogEventBegin(BV_OrthogonalizeVec,bv,0,0,0);
492: ksave = bv->k;
493: lsave = bv->l;
494: bv->l = -bv->nc; /* must also orthogonalize against constraints and leading columns */
495: if (!bv->buffer) { BVGetBufferVec(bv,&bv->buffer); }
496: BV_AllocateSignature(bv);
497: BVOrthogonalizeGS(bv,j,NULL,which,norm,lindep);
498: bv->k = ksave;
499: bv->l = lsave;
500: if (H) { BV_StoreCoefficients(bv,j,NULL,H); }
501: PetscLogEventEnd(BV_OrthogonalizeVec,bv,0,0,0);
502: PetscObjectStateIncrease((PetscObject)bv);
503: return(0);
504: }
506: /*
507: Block Gram-Schmidt: V2 = V2 - V1*R12, where R12 = V1'*V2
508: */
509: static PetscErrorCode BVOrthogonalize_BlockGS(BV V,Mat R)510: {
512: BV V1;
515: BVGetSplit(V,&V1,NULL);
516: BVDot(V,V1,R);
517: BVMult(V,-1.0,1.0,V1,R);
518: BVRestoreSplit(V,&V1,NULL);
519: return(0);
520: }
522: /*
523: Orthogonalize a set of vectors with Gram-Schmidt, column by column.
524: */
525: static PetscErrorCode BVOrthogonalize_GS(BV V,Mat R)526: {
528: PetscScalar *r=NULL;
529: PetscReal norm;
530: PetscInt j,ldr,lsave;
531: Vec v,w;
534: if (R) {
535: MatGetSize(R,&ldr,NULL);
536: MatDenseGetArray(R,&r);
537: }
538: if (V->matrix) {
539: BVGetCachedBV(V,&V->cached);
540: BVSetActiveColumns(V->cached,V->l,V->k);
541: }
542: for (j=V->l;j<V->k;j++) {
543: if (V->matrix && V->orthog_type==BV_ORTHOG_MGS) { /* fill cached BV */
544: BVGetColumn(V->cached,j,&v);
545: BVGetColumn(V,j,&w);
546: MatMult(V->matrix,w,v);
547: BVRestoreColumn(V,j,&w);
548: BVRestoreColumn(V->cached,j,&v);
549: }
550: if (R) {
551: BVOrthogonalizeColumn(V,j,NULL,&norm,NULL);
552: lsave = V->l;
553: V->l = -V->nc;
554: BV_StoreCoefficients(V,j,NULL,r+j*ldr);
555: V->l = lsave;
556: r[j+j*ldr] = norm;
557: } else {
558: BVOrthogonalizeColumn(V,j,NULL,&norm,NULL);
559: }
560: if (!norm) SETERRQ(PetscObjectComm((PetscObject)V),1,"Breakdown in BVOrthogonalize due to a linearly dependent column");
561: if (V->matrix && V->orthog_type==BV_ORTHOG_CGS) { /* fill cached BV */
562: BVGetColumn(V->cached,j,&v);
563: VecCopy(V->Bx,v);
564: BVRestoreColumn(V->cached,j,&v);
565: }
566: BVScaleColumn(V,j,1.0/norm);
567: }
568: if (R) { MatDenseRestoreArray(R,&r); }
569: return(0);
570: }
572: /*
573: BV_GetBufferMat - Create auxiliary seqdense matrix that wraps the bv->buffer.
574: */
575: PETSC_STATIC_INLINE PetscErrorCode BV_GetBufferMat(BV bv)576: {
578: PetscInt ld;
579: PetscScalar *array;
582: if (!bv->Abuffer) {
583: if (!bv->buffer) { BVGetBufferVec(bv,&bv->buffer); }
584: ld = bv->m+bv->nc;
585: VecGetArray(bv->buffer,&array);
586: MatCreateSeqDense(PETSC_COMM_SELF,ld,bv->m,array,&bv->Abuffer);
587: VecRestoreArray(bv->buffer,&array);
588: PetscLogObjectParent((PetscObject)bv,(PetscObject)bv->Abuffer);
589: }
590: return(0);
591: }
593: /*
594: BV_StoreCoeffsBlock_Default - Copy the contents of the BV buffer to a dense Mat
595: provided by the caller. Only columns l:k-1 are copied, restricting to the upper
596: triangular part if tri=PETSC_TRUE.
597: */
598: PETSC_STATIC_INLINE PetscErrorCode BV_StoreCoeffsBlock_Default(BV bv,Mat R,PetscBool tri)599: {
600: PetscErrorCode ierr;
601: const PetscScalar *bb;
602: PetscScalar *rr;
603: PetscInt j,ldr,ldb;
606: MatGetSize(R,&ldr,NULL);
607: MatDenseGetArray(R,&rr);
608: ldb = bv->m+bv->nc;
609: VecGetArrayRead(bv->buffer,&bb);
610: for (j=bv->l;j<bv->k;j++) {
611: PetscMemcpy(rr+j*ldr,bb+j*ldb,((tri?(j+1):bv->k)+bv->nc)*sizeof(PetscScalar));
612: }
613: VecRestoreArrayRead(bv->buffer,&bb);
614: MatDenseRestoreArray(R,&rr);
615: return(0);
616: }
618: /*
619: Orthogonalize a set of vectors with Cholesky: R=chol(V'*V), Q=V*inv(R)
620: */
621: static PetscErrorCode BVOrthogonalize_Chol(BV V,Mat Rin)622: {
624: Mat R,S;
627: BV_GetBufferMat(V);
628: R = V->Abuffer;
629: if (Rin) S = Rin; /* use Rin as a workspace for S */
630: else S = R;
631: if (V->l) { BVOrthogonalize_BlockGS(V,R); }
632: BVDot(V,V,R);
633: BVMatCholInv_LAPACK_Private(V,R,S);
634: BVMultInPlace(V,S,V->l,V->k);
635: if (Rin) { BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE); }
636: return(0);
637: }
639: /*
640: Orthogonalize a set of vectors with the Tall-Skinny QR method
641: */
642: static PetscErrorCode BVOrthogonalize_TSQR(BV V,Mat Rin)643: {
645: PetscScalar *pv,*r=NULL;
646: PetscInt ldr;
647: Mat R;
650: BV_GetBufferMat(V);
651: R = V->Abuffer;
652: if (V->l) { BVOrthogonalize_BlockGS(V,R); }
653: MatGetSize(R,&ldr,NULL);
654: MatDenseGetArray(R,&r);
655: BVGetArray(V,&pv);
656: BVOrthogonalize_LAPACK_TSQR(V,V->n,V->k-V->l,pv+(V->nc+V->l)*V->n,r+V->l*ldr+V->l,ldr);
657: BVRestoreArray(V,&pv);
658: MatDenseRestoreArray(R,&r);
659: if (Rin) { BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE); }
660: return(0);
661: }
663: /*
664: Orthogonalize a set of vectors with TSQR, but computing R only, then doing Q=V*inv(R)
665: */
666: static PetscErrorCode BVOrthogonalize_TSQRCHOL(BV V,Mat Rin)667: {
669: PetscScalar *pv,*r=NULL;
670: PetscInt ldr;
671: Mat R,S;
674: BV_GetBufferMat(V);
675: R = V->Abuffer;
676: if (Rin) S = Rin; /* use Rin as a workspace for S */
677: else S = R;
678: if (V->l) { BVOrthogonalize_BlockGS(V,R); }
679: MatGetSize(R,&ldr,NULL);
680: MatDenseGetArray(R,&r);
681: BVGetArray(V,&pv);
682: BVOrthogonalize_LAPACK_TSQR_OnlyR(V,V->n,V->k-V->l,pv+(V->nc+V->l)*V->n,r+V->l*ldr+V->l,ldr);
683: BVRestoreArray(V,&pv);
684: MatDenseRestoreArray(R,&r);
685: BVMatTriInv_LAPACK_Private(V,R,S);
686: BVMultInPlace(V,S,V->l,V->k);
687: if (Rin) { BV_StoreCoeffsBlock_Default(V,Rin,PETSC_TRUE); }
688: return(0);
689: }
691: /*
692: Orthogonalize a set of vectors with SVQB
693: */
694: static PetscErrorCode BVOrthogonalize_SVQB(BV V,Mat Rin)695: {
697: Mat R,S;
700: BV_GetBufferMat(V);
701: R = V->Abuffer;
702: if (Rin) S = Rin; /* use Rin as a workspace for S */
703: else S = R;
704: if (V->l) { BVOrthogonalize_BlockGS(V,R); }
705: BVDot(V,V,R);
706: BVMatSVQB_LAPACK_Private(V,R,S);
707: BVMultInPlace(V,S,V->l,V->k);
708: if (Rin) { BV_StoreCoeffsBlock_Default(V,Rin,PETSC_FALSE); }
709: return(0);
710: }
712: /*@
713: BVOrthogonalize - Orthogonalize all columns (starting from the leading ones),
714: that is, compute the QR decomposition.
716: Collective on BV718: Input Parameters:
719: + V - basis vectors to be orthogonalized (or B-orthogonalized)
720: - R - a sequential dense matrix (or NULL)
722: Output Parameters:
723: + V - the modified basis vectors
724: - R - (if not NULL) the triangular factor of the QR decomposition
726: Notes:
727: On input, matrix R must be a square sequential dense Mat, with at least as many
728: rows and columns as the number of active columns of V. The output satisfies
729: V0 = V*R (where V0 represent the input V) and V'*V = I (or V'*B*V = I if an
730: inner product matrix B has been specified with BVSetMatrix()).
732: If V has leading columns, then they are not modified (are assumed to be already
733: orthonormal) and the leading columns of R are not referenced. Let the
734: decomposition be
735: .vb
736: [ V01 V02 ] = [ V1 V2 ] [ R11 R12 ]
737: [ 0 R22 ]
738: .ve
739: then V1 is left unchanged (equal to V01) as well as R11 (it should satisfy
740: V01 = V1*R11).
742: Can pass NULL if R is not required.
744: The method to be used for block orthogonalization can be set with
745: BVSetOrthogonalization(). If set to GS, the computation is done column by
746: column with successive calls to BVOrthogonalizeColumn(). Note that in the
747: SVQB method the R factor is not upper triangular.
749: If V is rank-deficient or very ill-conditioned, that is, one or more columns are
750: (almost) linearly dependent with respect to the rest, then the algorithm may
751: break down or result in larger numerical error. Linearly dependent columns are
752: essentially replaced by random directions, and the corresponding diagonal entry
753: in R is set to (nearly) zero.
755: Level: intermediate
757: .seealso: BVOrthogonalizeColumn(), BVOrthogonalizeVec(), BVSetMatrix(), BVSetActiveColumns(), BVSetOrthogonalization(), BVOrthogBlockType758: @*/
759: PetscErrorCode BVOrthogonalize(BV V,Mat R)760: {
762: PetscBool match;
763: PetscInt m,n;
768: BVCheckSizes(V,1);
769: if (R) {
772: PetscObjectTypeCompare((PetscObject)R,MATSEQDENSE,&match);
773: if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Mat argument must be of type seqdense");
774: MatGetSize(R,&m,&n);
775: if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument is not square, it has %D rows and %D columns",m,n);
776: if (n<V->k) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat size %D is smaller than the number of BV active columns %D",n,V->k);
777: }
778: if (V->nc) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Not implemented for BV with constraints, use BVOrthogonalizeColumn() instead");
780: PetscLogEventBegin(BV_Orthogonalize,V,R,0,0);
781: switch (V->orthog_block) {
782: case BV_ORTHOG_BLOCK_GS: /* proceed column by column with Gram-Schmidt */
783: BVOrthogonalize_GS(V,R);
784: break;
785: case BV_ORTHOG_BLOCK_CHOL:
786: BVOrthogonalize_Chol(V,R);
787: break;
788: case BV_ORTHOG_BLOCK_TSQR:
789: if (V->matrix) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
790: BVOrthogonalize_TSQR(V,R);
791: break;
792: case BV_ORTHOG_BLOCK_TSQRCHOL:
793: if (V->matrix) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Orthogonalization method not available for non-standard inner product");
794: BVOrthogonalize_TSQRCHOL(V,R);
795: break;
796: case BV_ORTHOG_BLOCK_SVQB:
797: BVOrthogonalize_SVQB(V,R);
798: break;
799: }
800: PetscLogEventEnd(BV_Orthogonalize,V,R,0,0);
801: PetscObjectStateIncrease((PetscObject)V);
802: return(0);
803: }