Actual source code: vecutil.c

      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(MatMult(B,W[i],w));
     78:       else PetscCall(MatMult(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  - Hermitian 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^*BV$,
    129:    where $B$ is assumed to be Hermitian.

    131:    If `V`, `W` represent eigenvectors computed by SLEPc, this function will not work
    132:    correctly if one of the eigenvalues is complex when running with real scalars.

    134:    If `lev` is not `NULL`, it will contain the maximum entry of matrix
    135:    $W^*V - I$ (in absolute value) omitting the diagonal. Otherwise, the matrix $W^*V$
    136:    is printed.

    138:    Level: developer

    140: .seealso: `VecCheckOrthonormality()`
    141: @*/
    142: PetscErrorCode VecCheckOrthogonality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
    143: {
    144:   PetscFunctionBegin;
    145:   PetscAssertPointer(V,1);
    149:   if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
    150:   if (W) {
    151:     PetscAssertPointer(W,3);
    153:     PetscCheckSameComm(*V,1,*W,3);
    154:   }
    155:   PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_FALSE));
    156:   PetscFunctionReturn(PETSC_SUCCESS);
    157: }

    159: /*@
    160:    VecCheckOrthonormality - Checks (or prints) the level of (bi-)orthonormality
    161:    of a set of vectors.

    163:    Collective

    165:    Input Parameters:
    166: +  V  - a set of vectors
    167: .  nv - number of V vectors
    168: .  W  - an alternative set of vectors (optional)
    169: .  nw - number of W vectors
    170: .  B  - Hermitian matrix defining the inner product (optional)
    171: -  viewer - optional visualization context

    173:    Output Parameter:
    174: .  lev - level of orthogonality (optional)

    176:    Notes:
    177:    This function is equivalent to `VecCheckOrthonormality()`, but in addition it checks
    178:    that the diagonal of $W^*V$ (or $W^*BV$) is equal to all ones.

    180:    Level: developer

    182: .seealso: `VecCheckOrthogonality()`
    183: @*/
    184: PetscErrorCode VecCheckOrthonormality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
    185: {
    186:   PetscFunctionBegin;
    187:   PetscAssertPointer(V,1);
    191:   if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
    192:   if (W) {
    193:     PetscAssertPointer(W,3);
    195:     PetscCheckSameComm(*V,1,*W,3);
    196:   }
    197:   PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_TRUE));
    198:   PetscFunctionReturn(PETSC_SUCCESS);
    199: }

    201: /*@
    202:    VecDuplicateEmpty - Creates a new vector of the same type as an existing vector,
    203:    but without internal array.

    205:    Collective

    207:    Input Parameter:
    208: .  v - a vector to mimic

    210:    Output Parameter:
    211: .  newv - location to put new vector

    213:    Note:
    214:    This is similar to `VecDuplicate()`, but the new vector does not have an internal
    215:    array, so the intended usage is with `VecPlaceArray()`.

    217:    Level: developer

    219: .seealso: `VecDuplicate()`, `MatCreateVecsEmpty()`
    220: @*/
    221: PetscErrorCode VecDuplicateEmpty(Vec v,Vec *newv)
    222: {
    223:   PetscBool      standard,cuda,hip,mpi;
    224:   PetscInt       N,nloc,bs;

    226:   PetscFunctionBegin;
    228:   PetscAssertPointer(newv,2);

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

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

    260:    Logically Collective

    262:    Input Parameters:
    263: +  v    - the vector to be filled with random values
    264: .  rctx - the random number context (can be `NULL`)
    265: .  w1   - first work vector (can be `NULL`)
    266: -  w2   - second work vector (can be `NULL`)

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

    272:    Level: developer

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

    284:   PetscFunctionBegin;

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

    307:   PetscCall(VecSetRandom(w1,rctx));
    308:   PetscCall(VecSetRandom(w2,rctx));
    309:   PetscCall(VecGetLocalSize(v,&n));
    310:   PetscCall(VecGetArrayWrite(v,&z));
    311:   PetscCall(VecGetArrayRead(w1,&x));
    312:   PetscCall(VecGetArrayRead(w2,&y));
    313:   for (i=0;i<n;i++) {
    314: #if defined(PETSC_USE_COMPLEX)
    315:     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])));
    316: #else
    317:     z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
    318: #endif
    319:   }
    320:   PetscCall(VecRestoreArrayWrite(v,&z));
    321:   PetscCall(VecRestoreArrayRead(w1,&x));
    322:   PetscCall(VecRestoreArrayRead(w2,&y));

    324:   PetscCall(VecDestroy(&v1));
    325:   PetscCall(VecDestroy(&v2));
    326:   PetscCall(PetscRandomDestroy(&rand));
    327:   PetscFunctionReturn(PETSC_SUCCESS);
    328: }