Actual source code: vecutil.c

slepc-main 2024-11-09
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: */

 11: #include <slepc/private/vecimplslepc.h>

 13: /*@
 14:    VecNormalizeComplex - Normalizes a possibly complex vector by the 2-norm.

 16:    Collective

 18:    Input Parameters:
 19: +  xr - the real part of the vector (overwritten on output)
 20: .  xi - the imaginary part of the vector (not referenced if iscomplex is false)
 21: -  iscomplex - a flag indicating if the vector is complex

 23:    Output Parameter:
 24: .  norm - the vector norm before normalization (can be set to NULL)

 26:    Level: developer

 28: .seealso: BVNormalize()
 29: @*/
 30: PetscErrorCode VecNormalizeComplex(Vec xr,Vec xi,PetscBool iscomplex,PetscReal *norm)
 31: {
 32: #if !defined(PETSC_USE_COMPLEX)
 33:   PetscReal      normr,normi,alpha;
 34: #endif

 36:   PetscFunctionBegin;
 38: #if !defined(PETSC_USE_COMPLEX)
 39:   if (iscomplex) {
 41:     PetscCall(VecNormBegin(xr,NORM_2,&normr));
 42:     PetscCall(VecNormBegin(xi,NORM_2,&normi));
 43:     PetscCall(VecNormEnd(xr,NORM_2,&normr));
 44:     PetscCall(VecNormEnd(xi,NORM_2,&normi));
 45:     alpha = SlepcAbsEigenvalue(normr,normi);
 46:     if (norm) *norm = alpha;
 47:     alpha = 1.0 / alpha;
 48:     PetscCall(VecScale(xr,alpha));
 49:     PetscCall(VecScale(xi,alpha));
 50:   } else
 51: #endif
 52:     PetscCall(VecNormalize(xr,norm));
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: static PetscErrorCode VecCheckOrthogonality_Private(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev,PetscBool norm)
 57: {
 58:   PetscInt       i,j;
 59:   PetscScalar    *vals;
 60:   PetscBool      isascii;
 61:   Vec            w;

 63:   PetscFunctionBegin;
 64:   if (!lev) {
 65:     if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)*V),&viewer));
 67:     PetscCheckSameComm(*V,1,viewer,6);
 68:     PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
 69:     if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);
 70:   }

 72:   PetscCall(PetscMalloc1(nv,&vals));
 73:   if (B) PetscCall(VecDuplicate(V[0],&w));
 74:   if (lev) *lev = 0.0;
 75:   for (i=0;i<nw;i++) {
 76:     if (B) {
 77:       if (W) PetscCall(MatMultTranspose(B,W[i],w));
 78:       else PetscCall(MatMultTranspose(B,V[i],w));
 79:     } else {
 80:       if (W) w = W[i];
 81:       else w = V[i];
 82:     }
 83:     PetscCall(VecMDot(w,nv,V,vals));
 84:     for (j=0;j<nv;j++) {
 85:       if (lev) {
 86:         if (i!=j) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]));
 87:         else if (norm) {
 88:           if (PetscRealPart(vals[j])<0.0) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]+PetscRealConstant(1.0)));  /* indefinite case */
 89:           else *lev = PetscMax(*lev,PetscAbsScalar(vals[j]-PetscRealConstant(1.0)));
 90:         }
 91:       } else {
 92: #if !defined(PETSC_USE_COMPLEX)
 93:         PetscCall(PetscViewerASCIIPrintf(viewer," %12g  ",(double)vals[j]));
 94: #else
 95:         PetscCall(PetscViewerASCIIPrintf(viewer," %12g%+12gi ",(double)PetscRealPart(vals[j]),(double)PetscImaginaryPart(vals[j])));
 96: #endif
 97:       }
 98:     }
 99:     if (!lev) PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
100:   }
101:   PetscCall(PetscFree(vals));
102:   if (B) PetscCall(VecDestroy(&w));
103:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

106: /*@
107:    VecCheckOrthogonality - Checks (or prints) the level of (bi-)orthogonality
108:    of a set of vectors.

110:    Collective

112:    Input Parameters:
113: +  V  - a set of vectors
114: .  nv - number of V vectors
115: .  W  - an alternative set of vectors (optional)
116: .  nw - number of W vectors
117: .  B  - matrix defining the inner product (optional)
118: -  viewer - optional visualization context

120:    Output Parameter:
121: .  lev - level of orthogonality (optional)

123:    Notes:
124:    This function computes W'*V and prints the result. It is intended to check
125:    the level of bi-orthogonality of the vectors in the two sets. If W is equal
126:    to NULL then V is used, thus checking the orthogonality of the V vectors.

128:    If matrix B is provided then the check uses the B-inner product, W'*B*V.

130:    If lev is not NULL, it will contain the maximum entry of matrix
131:    W'*V - I (in absolute value) omitting the diagonal. Otherwise, the matrix W'*V
132:    is printed.

134:    Level: developer

136: .seealso: VecCheckOrthonormality()
137: @*/
138: PetscErrorCode VecCheckOrthogonality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
139: {
140:   PetscFunctionBegin;
141:   PetscAssertPointer(V,1);
145:   if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
146:   if (W) {
147:     PetscAssertPointer(W,3);
149:     PetscCheckSameComm(*V,1,*W,3);
150:   }
151:   PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_FALSE));
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: /*@
156:    VecCheckOrthonormality - Checks (or prints) the level of (bi-)orthonormality
157:    of a set of vectors.

159:    Collective

161:    Input Parameters:
162: +  V  - a set of vectors
163: .  nv - number of V vectors
164: .  W  - an alternative set of vectors (optional)
165: .  nw - number of W vectors
166: .  B  - matrix defining the inner product (optional)
167: -  viewer - optional visualization context

169:    Output Parameter:
170: .  lev - level of orthogonality (optional)

172:    Notes:
173:    This function is equivalent to VecCheckOrthonormality(), but in addition it checks
174:    that the diagonal of W'*V (or W'*B*V) is equal to all ones.

176:    Level: developer

178: .seealso: VecCheckOrthogonality()
179: @*/
180: PetscErrorCode VecCheckOrthonormality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
181: {
182:   PetscFunctionBegin;
183:   PetscAssertPointer(V,1);
187:   if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
188:   if (W) {
189:     PetscAssertPointer(W,3);
191:     PetscCheckSameComm(*V,1,*W,3);
192:   }
193:   PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_TRUE));
194:   PetscFunctionReturn(PETSC_SUCCESS);
195: }

197: /*@
198:    VecDuplicateEmpty - Creates a new vector of the same type as an existing vector,
199:    but without internal array.

201:    Collective

203:    Input Parameters:
204: .  v - a vector to mimic

206:    Output Parameter:
207: .  newv - location to put new vector

209:    Note:
210:    This is similar to VecDuplicate(), but the new vector does not have an internal
211:    array, so the intended usage is with VecPlaceArray().

213:    Level: developer

215: .seealso: MatCreateVecsEmpty()
216: @*/
217: PetscErrorCode VecDuplicateEmpty(Vec v,Vec *newv)
218: {
219:   PetscBool      standard,cuda,hip,mpi;
220:   PetscInt       N,nloc,bs;

222:   PetscFunctionBegin;
224:   PetscAssertPointer(newv,2);

227:   PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&standard,VECSEQ,VECMPI,""));
228:   PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,""));
229:   PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&hip,VECSEQHIP,VECMPIHIP,""));
230:   if (standard || cuda || hip) {
231:     PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
232:     PetscCall(VecGetLocalSize(v,&nloc));
233:     PetscCall(VecGetSize(v,&N));
234:     PetscCall(VecGetBlockSize(v,&bs));
235:     if (cuda) {
236: #if defined(PETSC_HAVE_CUDA)
237:       if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
238:       else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
239: #endif
240:     } else if (hip) {
241: #if defined(PETSC_HAVE_HIP)
242:       if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
243:       else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
244: #endif
245:     } else {
246:       if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
247:       else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
248:     }
249:   } else PetscCall(VecDuplicate(v,newv)); /* standard duplicate, with internal array */
250:   PetscFunctionReturn(PETSC_SUCCESS);
251: }

253: /*@
254:    VecSetRandomNormal - Sets all components of a vector to normally distributed random values.

256:    Logically Collective

258:    Input Parameters:
259: +  v    - the vector to be filled with random values
260: .  rctx - the random number context (can be NULL)
261: .  w1   - first work vector (can be NULL)
262: -  w2   - second work vector (can be NULL)

264:    Notes:
265:    Fills the two work vectors with uniformly distributed random values (VecSetRandom)
266:    and then applies the Box-Muller transform to get normally distributed values on v.

268:    Level: developer

270: .seealso: VecSetRandom()
271: @*/
272: PetscErrorCode VecSetRandomNormal(Vec v,PetscRandom rctx,Vec w1,Vec w2)
273: {
274:   const PetscScalar *x,*y;
275:   PetscScalar       *z;
276:   PetscInt          n,i;
277:   PetscRandom       rand=NULL;
278:   Vec               v1=NULL,v2=NULL;

280:   PetscFunctionBegin;

287:   if (!rctx) {
288:     PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)v),&rand));
289:     PetscCall(PetscRandomSetFromOptions(rand));
290:     rctx = rand;
291:   }
292:   if (!w1) {
293:     PetscCall(VecDuplicate(v,&v1));
294:     w1 = v1;
295:   }
296:   if (!w2) {
297:     PetscCall(VecDuplicate(v,&v2));
298:     w2 = v2;
299:   }
300:   PetscCheckSameTypeAndComm(v,1,w1,3);
301:   PetscCheckSameTypeAndComm(v,1,w2,4);

303:   PetscCall(VecSetRandom(w1,rctx));
304:   PetscCall(VecSetRandom(w2,rctx));
305:   PetscCall(VecGetLocalSize(v,&n));
306:   PetscCall(VecGetArrayWrite(v,&z));
307:   PetscCall(VecGetArrayRead(w1,&x));
308:   PetscCall(VecGetArrayRead(w2,&y));
309:   for (i=0;i<n;i++) {
310: #if defined(PETSC_USE_COMPLEX)
311:     z[i] = PetscCMPLX(PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(y[i])),PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(y[i])));
312: #else
313:     z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
314: #endif
315:   }
316:   PetscCall(VecRestoreArrayWrite(v,&z));
317:   PetscCall(VecRestoreArrayRead(w1,&x));
318:   PetscCall(VecRestoreArrayRead(w2,&y));

320:   PetscCall(VecDestroy(&v1));
321:   PetscCall(VecDestroy(&v2));
322:   if (!rctx) PetscCall(PetscRandomDestroy(&rand));
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }