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: }