Actual source code: vecutil.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: */

 11: #include <slepc/private/vecimplslepc.h>       /*I "slepcvec.h" I*/

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

 16:    Collective on Vec

 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
 27: @*/
 28: PetscErrorCode VecNormalizeComplex(Vec xr,Vec xi,PetscBool iscomplex,PetscReal *norm)
 29: {
 31: #if !defined(PETSC_USE_COMPLEX)
 32:   PetscReal      normr,normi,alpha;
 33: #endif

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

 57: /*@C
 58:    VecCheckOrthogonality - Checks (or prints) the level of (bi-)orthogonality
 59:    of a set of vectors.

 61:    Collective on Vec

 63:    Input parameters:
 64: +  V  - a set of vectors
 65: .  nv - number of V vectors
 66: .  W  - an alternative set of vectors (optional)
 67: .  nw - number of W vectors
 68: .  B  - matrix defining the inner product (optional)
 69: -  viewer - optional visualization context

 71:    Output parameter:
 72: .  lev - level of orthogonality (optional)

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

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

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

 85:    Level: developer
 86: @*/
 87: PetscErrorCode VecCheckOrthogonality(Vec *V,PetscInt nv,Vec *W,PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
 88: {
 90:   PetscInt       i,j;
 91:   PetscScalar    *vals;
 92:   PetscBool      isascii;
 93:   Vec            w;

 96:   if (nv<=0 || nw<=0) return(0);
 99:   if (!lev) {
100:     if (!viewer) {
101:       PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)*V),&viewer);
102:     }
105:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
106:     if (!isascii) return(0);
107:   }

109:   PetscMalloc1(nv,&vals);
110:   if (B) {
111:     VecDuplicate(V[0],&w);
112:   }
113:   if (lev) *lev = 0.0;
114:   for (i=0;i<nw;i++) {
115:     if (B) {
116:       if (W) {
117:         MatMultTranspose(B,W[i],w);
118:       } else {
119:         MatMultTranspose(B,V[i],w);
120:       }
121:     } else {
122:       if (W) w = W[i];
123:       else w = V[i];
124:     }
125:     VecMDot(w,nv,V,vals);
126:     for (j=0;j<nv;j++) {
127:       if (lev) {
128:         if (i!=j) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]));
129:       } else {
130: #if !defined(PETSC_USE_COMPLEX)
131:         PetscViewerASCIIPrintf(viewer," %12g  ",(double)vals[j]);
132: #else
133:         PetscViewerASCIIPrintf(viewer," %12g%+12gi ",(double)PetscRealPart(vals[j]),(double)PetscImaginaryPart(vals[j]));
134: #endif
135:       }
136:     }
137:     if (!lev) { PetscViewerASCIIPrintf(viewer,"\n"); }
138:   }
139:   PetscFree(vals);
140:   if (B) {
141:     VecDestroy(&w);
142:   }
143:   return(0);
144: }

146: /*@
147:    VecDuplicateEmpty - Creates a new vector of the same type as an existing vector,
148:    but without internal array.

150:    Collective on Vec

152:    Input parameters:
153: .  v - a vector to mimic

155:    Output parameter:
156: .  newv - location to put new vector

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

162:    Level: developer
163: @*/
164: PetscErrorCode VecDuplicateEmpty(Vec v,Vec *newv)
165: {
167:   PetscBool      sup,cuda,mpi;
168:   PetscInt       N,nloc,bs;


175:   PetscObjectTypeCompareAny((PetscObject)v,&sup,VECSEQ,VECMPI,VECSEQCUDA,VECMPICUDA,"");
176:   if (!sup) SETERRQ1(PetscObjectComm((PetscObject)v),PETSC_ERR_SUP,"Vector type %s not supported",((PetscObject)v)->type_name);
177:   PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,"");
178:   PetscObjectTypeCompareAny((PetscObject)v,&mpi,VECMPI,VECMPICUDA,"");
179:   VecGetLocalSize(v,&nloc);
180:   VecGetSize(v,&N);
181:   VecGetBlockSize(v,&bs);

183:   if (cuda) {
184: #if defined(PETSC_HAVE_CUDA)
185:     if (mpi) {
186:       VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv);
187:     } else {
188:       VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv);
189:     }
190: #endif
191:   } else {
192:     if (mpi) {
193:       VecCreateMPIWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv);
194:     } else {
195:       VecCreateSeqWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv);
196:     }
197:   }
198:   return(0);
199: }