Actual source code: bvbiorthog.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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 bi-orthogonalization routines
 12: */

 14: #include <slepc/private/bvimpl.h>          /*I   "slepcbv.h"   I*/

 16: /*
 17:    BVBiorthogonalizeMGS1 - Compute one step of Modified Gram-Schmidt bi-orthogonalization
 18: */
 19: static PetscErrorCode BVBiorthogonalizeMGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
 20: {
 22:   PetscInt       i;
 23:   PetscScalar    dot;
 24:   Vec            vi,wi;

 27:   for (i=-V->nc;i<V->k;i++) {
 28:     BVGetColumn(W,i,&wi);
 29:     /* h_i = ( v, w_i ) */
 30:     VecDot(v,wi,&dot);
 31:     BVRestoreColumn(W,i,&wi);
 32:     /* v <- v - h_i v_i */
 33:     BV_SetValue(V,i,0,c,dot);
 34:     BVGetColumn(V,i,&vi);
 35:     VecAXPY(v,-dot,vi);
 36:     BVRestoreColumn(V,i,&vi);
 37:   }
 38:   BV_AddCoefficients(V,V->k,h,c);
 39:   return(0);
 40: }

 42: /*
 43:    BVBiorthogonalizeCGS1 - Compute one step of CGS bi-orthogonalization: v = (I-V*W')*v
 44: */
 45: static PetscErrorCode BVBiorthogonalizeCGS1(BV V,BV W,Vec v,PetscScalar *h,PetscScalar *c)
 46: {

 50:   /* h = W'*v */
 51:   BVDotVec(W,v,c);

 53:   /* v = v - V h */
 54:   BVMultVec(V,-1.0,1.0,v,c);

 56:   BV_AddCoefficients(V,V->k,h,c);
 57:   return(0);
 58: }

 60: #define BVBiorthogonalizeGS1(a,b,c,d,e) ((V->orthog_type==BV_ORTHOG_MGS)?BVBiorthogonalizeMGS1:BVBiorthogonalizeCGS1)(a,b,c,d,e)

 62: /*
 63:    BVBiorthogonalizeGS - Orthogonalize with (classical or modified) Gram-Schmidt

 65:    V, W - the two basis vectors objects
 66:    v    - the vector to bi-orthogonalize
 67: */
 68: static PetscErrorCode BVBiorthogonalizeGS(BV V,BV W,Vec v)
 69: {
 71:   PetscScalar    *h,*c;

 74:   h = V->h;
 75:   c = V->c;
 76:   BV_CleanCoefficients(V,V->k,h);
 77:   BVBiorthogonalizeGS1(V,W,v,h,c);
 78:   if (V->orthog_ref!=BV_ORTHOG_REFINE_NEVER) {
 79:     BVBiorthogonalizeGS1(V,W,v,h,c);
 80:   }
 81:   return(0);
 82: }

 84: /*@
 85:    BVBiorthogonalizeColumn - Bi-orthogonalize a column of two BV objects.

 87:    Collective on BV

 89:    Input Parameters:
 90: +  V,W - two basis vectors contexts
 91: -  j   - index of column to be bi-orthonormalized

 93:    Notes:
 94:    This function bi-orthogonalizes vectors V[j],W[j] against W[0..j-1],
 95:    and V[0..j-1], respectively, so that W[0..j]'*V[0..j] = diagonal.

 97:    Level: advanced

 99: .seealso: BVOrthogonalizeColumn(), BVBiorthonormalizeColumn()
100: @*/
101: PetscErrorCode BVBiorthogonalizeColumn(BV V,BV W,PetscInt j)
102: {
104:   PetscInt       ksavev,lsavev,ksavew,lsavew;
105:   Vec            y,z;

112:   BVCheckSizes(V,1);
114:   BVCheckSizes(W,2);
116:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
117:   if (j>=V->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but V only has %D columns",j,V->m);
118:   if (j>=W->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but W only has %D columns",j,W->m);
119:   if (V->n!=W->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, W %D",V->n,W->n);
120:   if (V->matrix || W->matrix) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
121:   if (V->nc || W->nc) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
122:   if (V->ops->gramschmidt || W->ops->gramschmidt) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");

124:   /* bi-orthogonalize */
125:   PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0);
126:   ksavev = V->k;
127:   lsavev = V->l;
128:   ksavew = W->k;
129:   lsavew = W->l;
130:   V->k = j;
131:   V->l = -V->nc;  /* must also bi-orthogonalize against constraints and leading columns */
132:   W->k = j;
133:   W->l = -W->nc;
134:   BV_AllocateCoeffs(V);
135:   BV_AllocateCoeffs(W);
136:   BVGetColumn(V,j,&y);
137:   BVBiorthogonalizeGS(V,W,y);
138:   BVRestoreColumn(V,j,&y);
139:   BVGetColumn(W,j,&z);
140:   BVBiorthogonalizeGS(W,V,z);
141:   BVRestoreColumn(W,j,&z);
142:   V->k = ksavev;
143:   V->l = lsavev;
144:   W->k = ksavew;
145:   W->l = lsavew;
146:   PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0);
147:   PetscObjectStateIncrease((PetscObject)V);
148:   PetscObjectStateIncrease((PetscObject)W);
149:   return(0);
150: }
151: /*@
152:    BVBiorthonormalizeColumn - Bi-orthonormalize a column of two BV objects.

154:    Collective on BV

156:    Input Parameters:
157: +  V,W - two basis vectors contexts
158: -  j   - index of column to be bi-orthonormalized

160:    Output Parameters:
161: .  delta - (optional) value used for normalization

163:    Notes:
164:    This function first bi-orthogonalizes vectors V[j],W[j] against W[0..j-1],
165:    and V[0..j-1], respectively. Then, it scales the vectors with 1/delta, so
166:    that the resulting vectors satisfy W[j]'*V[j] = 1.

168:    Level: advanced

170: .seealso: BVOrthonormalizeColumn(), BVBiorthogonalizeColumn()
171: @*/
172: PetscErrorCode BVBiorthonormalizeColumn(BV V,BV W,PetscInt j,PetscReal *delta)
173: {
175:   PetscScalar    alpha;
176:   PetscReal      deltat;
177:   PetscInt       ksavev,lsavev,ksavew,lsavew;
178:   Vec            y,z;

185:   BVCheckSizes(V,1);
187:   BVCheckSizes(W,2);
189:   if (j<0) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
190:   if (j>=V->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but V only has %D columns",j,V->m);
191:   if (j>=W->m) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%D but W only has %D columns",j,W->m);
192:   if (V->n!=W->n) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %D, W %D",V->n,W->n);
193:   if (V->matrix || W->matrix) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W must not have an inner product matrix");
194:   if (V->nc || W->nc) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_WRONGSTATE,"V,W cannot have different number of constraints");
195:   if (V->ops->gramschmidt || W->ops->gramschmidt) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Object has a special GS function");

197:   /* bi-orthogonalize */
198:   PetscLogEventBegin(BV_OrthogonalizeVec,V,0,0,0);
199:   ksavev = V->k;
200:   lsavev = V->l;
201:   ksavew = W->k;
202:   lsavew = W->l;
203:   V->k = j;
204:   V->l = -V->nc;  /* must also bi-orthogonalize against constraints and leading columns */
205:   W->k = j;
206:   W->l = -W->nc;
207:   BV_AllocateCoeffs(V);
208:   BV_AllocateCoeffs(W);
209:   BVGetColumn(V,j,&y);
210:   BVBiorthogonalizeGS(V,W,y);
211:   BVRestoreColumn(V,j,&y);
212:   BVGetColumn(W,j,&z);
213:   BVBiorthogonalizeGS(W,V,z);
214:   BVRestoreColumn(W,j,&z);
215:   V->k = ksavev;
216:   V->l = lsavev;
217:   W->k = ksavew;
218:   W->l = lsavew;
219:   PetscLogEventEnd(BV_OrthogonalizeVec,V,0,0,0);

221:   /* scale */
222:   PetscLogEventBegin(BV_Scale,V,0,0,0);
223:   BVGetColumn(V,j,&y);
224:   BVGetColumn(W,j,&z);
225:   VecDot(y,z,&alpha);
226:   BVRestoreColumn(V,j,&y);
227:   BVRestoreColumn(W,j,&z);
228:   deltat = PetscSqrtReal(PetscAbsScalar(alpha));
229:   if (V->n) { (*V->ops->scale)(V,j,1.0/PetscConj(alpha/deltat)); }
230:   if (W->n) { (*W->ops->scale)(W,j,1.0/deltat); }
231:   PetscLogEventEnd(BV_Scale,V,0,0,0);
232:   if (delta) *delta = deltat;
233:   PetscObjectStateIncrease((PetscObject)V);
234:   PetscObjectStateIncrease((PetscObject)W);
235:   return(0);
236: }