Actual source code: test11.c
slepc-3.11.2 2019-07-30
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: static char help[] = "Test BV block orthogonalization.\n\n";
13: #include <slepcbv.h>
15: /*
16: Compute the Frobenius norm ||A(l:k,l:k)-diag||_F
17: */
18: PetscErrorCode MyMatNorm(Mat A,PetscInt lda,PetscInt l,PetscInt k,PetscScalar diag,PetscReal *norm)
19: {
21: PetscInt i,j;
22: PetscScalar *pA;
23: PetscReal s,val;
26: MatDenseGetArray(A,&pA);
27: s = 0.0;
28: for (i=l;i<k;i++) {
29: for (j=l;j<k;j++) {
30: val = (i==j)? PetscAbsScalar(pA[i+j*lda]-diag): PetscAbsScalar(pA[i+j*lda]);
31: s += val*val;
32: }
33: }
34: *norm = PetscSqrtReal(s);
35: MatDenseRestoreArray(A,&pA);
36: return(0);
37: }
39: int main(int argc,char **argv)
40: {
42: BV X,Y,Z,cached;
43: Mat B=NULL,M,R=NULL;
44: Vec v,t;
45: PetscInt i,j,n=20,l=2,k=8,Istart,Iend;
46: PetscViewer view;
47: PetscBool withb,resid,rand,verbose;
48: PetscReal norm;
49: PetscScalar alpha;
51: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
52: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
53: PetscOptionsGetInt(NULL,NULL,"-l",&l,NULL);
54: PetscOptionsGetInt(NULL,NULL,"-k",&k,NULL);
55: PetscOptionsHasName(NULL,NULL,"-withb",&withb);
56: PetscOptionsHasName(NULL,NULL,"-resid",&resid);
57: PetscOptionsHasName(NULL,NULL,"-rand",&rand);
58: PetscOptionsHasName(NULL,NULL,"-verbose",&verbose);
59: PetscPrintf(PETSC_COMM_WORLD,"Test BV block orthogonalization (length %D, l=%D, k=%D)%s.\n",n,l,k,withb?" with non-standard inner product":"");
61: /* Create template vector */
62: VecCreate(PETSC_COMM_WORLD,&t);
63: VecSetSizes(t,PETSC_DECIDE,n);
64: VecSetFromOptions(t);
66: /* Create BV object X */
67: BVCreate(PETSC_COMM_WORLD,&X);
68: PetscObjectSetName((PetscObject)X,"X");
69: BVSetSizesFromVec(X,t,k);
70: BVSetFromOptions(X);
72: /* Set up viewer */
73: PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&view);
74: if (verbose) {
75: PetscViewerPushFormat(view,PETSC_VIEWER_ASCII_MATLAB);
76: }
78: /* Fill X entries */
79: if (rand) {
80: BVSetRandom(X);
81: } else {
82: for (j=0;j<k;j++) {
83: BVGetColumn(X,j,&v);
84: VecSet(v,0.0);
85: for (i=0;i<=n/2;i++) {
86: if (i+j<n) {
87: alpha = (3.0*i+j-2)/(2*(i+j+1));
88: VecSetValue(v,i+j,alpha,INSERT_VALUES);
89: }
90: }
91: VecAssemblyBegin(v);
92: VecAssemblyEnd(v);
93: BVRestoreColumn(X,j,&v);
94: }
95: }
96: if (verbose) {
97: BVView(X,view);
98: }
100: if (withb) {
101: /* Create inner product matrix */
102: MatCreate(PETSC_COMM_WORLD,&B);
103: MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n);
104: MatSetFromOptions(B);
105: MatSetUp(B);
106: PetscObjectSetName((PetscObject)B,"B");
108: MatGetOwnershipRange(B,&Istart,&Iend);
109: for (i=Istart;i<Iend;i++) {
110: if (i>0) { MatSetValue(B,i,i-1,-1.0,INSERT_VALUES); }
111: if (i<n-1) { MatSetValue(B,i,i+1,-1.0,INSERT_VALUES); }
112: MatSetValue(B,i,i,2.0,INSERT_VALUES);
113: }
114: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
115: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
116: if (verbose) {
117: MatView(B,view);
118: }
119: BVSetMatrix(X,B,PETSC_FALSE);
120: }
122: /* Create copy on Y */
123: BVDuplicate(X,&Y);
124: PetscObjectSetName((PetscObject)Y,"Y");
125: BVCopy(X,Y);
126: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&M);
128: if (resid) {
129: /* Create matrix R to store triangular factor */
130: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&R);
131: PetscObjectSetName((PetscObject)R,"R");
132: }
134: if (l>0) {
135: /* First orthogonalize leading columns */
136: PetscPrintf(PETSC_COMM_WORLD,"Orthogonalizing leading columns\n");
137: BVSetActiveColumns(Y,0,l);
138: BVSetActiveColumns(X,0,l);
139: BVOrthogonalize(Y,R);
140: if (verbose) {
141: BVView(Y,view);
142: if (resid) { MatView(R,view); }
143: }
145: if (withb) {
146: /* Extract cached BV and check it is equal to B*X */
147: BVGetCachedBV(Y,&cached);
148: BVDuplicate(X,&Z);
149: BVSetMatrix(Z,NULL,PETSC_FALSE);
150: BVSetActiveColumns(Z,0,l);
151: BVCopy(X,Z);
152: BVMatMult(X,B,Z);
153: BVMult(Z,-1.0,1.0,cached,NULL);
154: BVNorm(Z,NORM_FROBENIUS,&norm);
155: if (norm<100*PETSC_MACHINE_EPSILON) {
156: PetscPrintf(PETSC_COMM_WORLD," Difference ||cached-BX|| < 100*eps\n");
157: } else {
158: PetscPrintf(PETSC_COMM_WORLD," Difference ||cached-BX||: %g\n",(double)norm);
159: }
160: BVDestroy(&Z);
161: }
163: /* Check orthogonality */
164: BVDot(Y,Y,M);
165: MyMatNorm(M,k,0,l,1.0,&norm);
166: if (norm<100*PETSC_MACHINE_EPSILON) {
167: PetscPrintf(PETSC_COMM_WORLD," Level of orthogonality of Q1 < 100*eps\n");
168: } else {
169: PetscPrintf(PETSC_COMM_WORLD," Level of orthogonality of Q1: %g\n",(double)norm);
170: }
172: if (resid) {
173: /* Check residual */
174: BVDuplicate(X,&Z);
175: BVSetMatrix(Z,NULL,PETSC_FALSE);
176: BVSetActiveColumns(Z,0,l);
177: BVCopy(X,Z);
178: BVMult(Z,-1.0,1.0,Y,R);
179: BVNorm(Z,NORM_FROBENIUS,&norm);
180: if (norm<100*PETSC_MACHINE_EPSILON) {
181: PetscPrintf(PETSC_COMM_WORLD," Residual ||X1-Q1*R11|| < 100*eps\n");
182: } else {
183: PetscPrintf(PETSC_COMM_WORLD," Residual ||X1-Q1*R11||: %g\n",(double)norm);
184: }
185: BVDestroy(&Z);
186: }
188: }
190: /* Now orthogonalize the rest of columns */
191: PetscPrintf(PETSC_COMM_WORLD,"Orthogonalizing active columns\n");
192: BVSetActiveColumns(Y,l,k);
193: BVSetActiveColumns(X,l,k);
194: BVOrthogonalize(Y,R);
195: if (verbose) {
196: BVView(Y,view);
197: if (resid) { MatView(R,view); }
198: }
200: if (l>0) {
201: /* Check orthogonality */
202: BVDot(Y,Y,M);
203: MyMatNorm(M,k,l,k,1.0,&norm);
204: if (norm<100*PETSC_MACHINE_EPSILON) {
205: PetscPrintf(PETSC_COMM_WORLD," Level of orthogonality of Q2 < 100*eps\n");
206: } else {
207: PetscPrintf(PETSC_COMM_WORLD," Level of orthogonality of Q2: %g\n",(double)norm);
208: }
209: }
211: /* Check the complete decomposition */
212: PetscPrintf(PETSC_COMM_WORLD,"Overall decomposition\n");
213: BVSetActiveColumns(Y,0,k);
214: BVSetActiveColumns(X,0,k);
216: /* Check orthogonality */
217: BVDot(Y,Y,M);
218: MyMatNorm(M,k,0,k,1.0,&norm);
219: if (norm<100*PETSC_MACHINE_EPSILON) {
220: PetscPrintf(PETSC_COMM_WORLD," Level of orthogonality of Q < 100*eps\n");
221: } else {
222: PetscPrintf(PETSC_COMM_WORLD," Level of orthogonality of Q: %g\n",(double)norm);
223: }
225: if (resid) {
226: /* Check residual */
227: BVMult(X,-1.0,1.0,Y,R);
228: BVSetMatrix(X,NULL,PETSC_FALSE);
229: BVNorm(X,NORM_FROBENIUS,&norm);
230: if (norm<100*PETSC_MACHINE_EPSILON) {
231: PetscPrintf(PETSC_COMM_WORLD," Residual ||X-Q*R|| < 100*eps\n");
232: } else {
233: PetscPrintf(PETSC_COMM_WORLD," Residual ||X-Q*R||: %g\n",(double)norm);
234: }
235: MatDestroy(&R);
236: }
238: if (B) { MatDestroy(&B); }
239: MatDestroy(&M);
240: BVDestroy(&X);
241: BVDestroy(&Y);
242: VecDestroy(&t);
243: SlepcFinalize();
244: return ierr;
245: }
247: /*TEST
249: test:
250: suffix: 1
251: nsize: 2
252: args: -bv_orthog_block gs -bv_type {{vecs contiguous svec mat}shared output}
253: output_file: output/test11_1.out
255: test:
256: suffix: 1_cuda
257: nsize: 2
258: args: -bv_orthog_block gs -bv_type svec -vec_type cuda
259: requires: cuda
260: output_file: output/test11_1.out
262: test:
263: suffix: 2
264: nsize: 2
265: args: -bv_orthog_block chol -bv_type {{vecs contiguous svec mat}shared output}
266: output_file: output/test11_1.out
268: test:
269: suffix: 2_cuda
270: nsize: 2
271: args: -bv_orthog_block chol -bv_type svec -vec_type cuda
272: requires: cuda
273: output_file: output/test11_1.out
275: test:
276: suffix: 3
277: nsize: 2
278: args: -bv_orthog_block tsqr -bv_type {{vecs contiguous svec mat}shared output}
279: output_file: output/test11_1.out
281: test:
282: suffix: 3_cuda
283: nsize: 2
284: args: -bv_orthog_block tsqr -bv_type svec -vec_type cuda
285: requires: cuda
286: output_file: output/test11_1.out
288: test:
289: suffix: 4
290: nsize: 2
291: args: -withb -bv_orthog_block gs -bv_type {{vecs contiguous svec mat}shared output}
292: output_file: output/test11_4.out
294: test:
295: suffix: 4_cuda
296: nsize: 2
297: args: -withb -bv_orthog_block gs -bv_type svec -vec_type cuda -mat_type aijcusparse
298: requires: cuda
299: output_file: output/test11_4.out
301: test:
302: suffix: 5
303: nsize: 2
304: args: -withb -bv_orthog_block chol -bv_type {{vecs contiguous svec mat}shared output}
305: output_file: output/test11_4.out
307: test:
308: suffix: 5_cuda
309: nsize: 2
310: args: -withb -bv_orthog_block chol -bv_type svec -vec_type cuda -mat_type aijcusparse
311: requires: cuda
312: output_file: output/test11_4.out
314: test:
315: suffix: 6
316: nsize: 2
317: args: -resid -bv_orthog_block gs -bv_type {{vecs contiguous svec mat}shared output}
318: output_file: output/test11_6.out
320: test:
321: suffix: 6_cuda
322: nsize: 2
323: args: -resid -bv_orthog_block gs -bv_type svec -vec_type cuda
324: requires: cuda
325: output_file: output/test11_6.out
327: test:
328: suffix: 7
329: nsize: 2
330: args: -resid -bv_orthog_block chol -bv_type {{vecs contiguous svec mat}shared output}
331: output_file: output/test11_6.out
333: test:
334: suffix: 7_cuda
335: nsize: 2
336: args: -resid -bv_orthog_block chol -bv_type svec -vec_type cuda
337: requires: cuda
338: output_file: output/test11_6.out
340: test:
341: suffix: 8
342: nsize: 2
343: args: -resid -bv_orthog_block tsqr -bv_type {{vecs contiguous svec mat}shared output}
344: output_file: output/test11_6.out
346: test:
347: suffix: 8_cuda
348: nsize: 2
349: args: -resid -bv_orthog_block tsqr -bv_type svec -vec_type cuda
350: requires: cuda
351: output_file: output/test11_6.out
353: test:
354: suffix: 9
355: nsize: 2
356: args: -resid -withb -bv_orthog_block gs -bv_type {{vecs contiguous svec mat}shared output}
357: output_file: output/test11_9.out
359: test:
360: suffix: 9_cuda
361: nsize: 2
362: args: -resid -withb -bv_orthog_block gs -bv_type svec -vec_type cuda -mat_type aijcusparse
363: requires: cuda
364: output_file: output/test11_9.out
366: test:
367: suffix: 10
368: nsize: 2
369: args: -resid -withb -bv_orthog_block chol -bv_type {{vecs contiguous svec mat}shared output}
370: output_file: output/test11_9.out
372: test:
373: suffix: 10_cuda
374: nsize: 2
375: args: -resid -withb -bv_orthog_block chol -bv_type svec -vec_type cuda -mat_type aijcusparse
376: requires: cuda
377: output_file: output/test11_9.out
379: test:
380: suffix: 11
381: nsize: 7
382: args: -bv_orthog_block tsqr -bv_type {{vecs contiguous svec mat}shared output}
383: output_file: output/test11_1.out
385: test:
386: suffix: 11_cuda
387: nsize: 7
388: args: -bv_orthog_block tsqr -bv_type svec -vec_type cuda
389: requires: cuda
390: output_file: output/test11_1.out
392: test:
393: suffix: 12
394: nsize: 9
395: args: -resid -n 180 -l 0 -k 7 -bv_orthog_block tsqr -bv_type {{vecs contiguous svec mat}shared output}
396: requires: !single
397: output_file: output/test11_12.out
399: test:
400: suffix: 12_cuda
401: nsize: 9
402: args: -resid -n 180 -l 0 -k 7 -bv_orthog_block tsqr -bv_type svec -vec_type cuda
403: requires: !single cuda
404: output_file: output/test11_12.out
406: test:
407: suffix: 13
408: nsize: 2
409: args: -bv_orthog_block tsqrchol -bv_type {{vecs contiguous svec mat}shared output}
410: output_file: output/test11_1.out
412: test:
413: suffix: 13_cuda
414: nsize: 2
415: args: -bv_orthog_block tsqrchol -bv_type svec -vec_type cuda
416: requires: cuda
417: output_file: output/test11_1.out
419: test:
420: suffix: 14
421: nsize: 2
422: args: -resid -bv_orthog_block tsqrchol -bv_type {{vecs contiguous svec mat}shared output}
423: output_file: output/test11_6.out
425: test:
426: suffix: 14_cuda
427: nsize: 2
428: args: -resid -bv_orthog_block tsqrchol -bv_type svec -vec_type cuda
429: requires: cuda
430: output_file: output/test11_6.out
432: test:
433: suffix: 15
434: nsize: 2
435: args: -bv_orthog_block svqb -bv_type {{vecs contiguous svec mat}shared output}
436: output_file: output/test11_1.out
438: test:
439: suffix: 15_cuda
440: nsize: 2
441: args: -bv_orthog_block svqb -bv_type svec -vec_type cuda
442: requires: cuda
443: output_file: output/test11_1.out
445: test:
446: suffix: 16
447: nsize: 2
448: args: -withb -bv_orthog_block svqb -bv_type {{vecs contiguous svec mat}shared output}
449: output_file: output/test11_4.out
451: test:
452: suffix: 16_cuda
453: nsize: 2
454: args: -withb -bv_orthog_block svqb -bv_type svec -vec_type cuda
455: requires: cuda
456: output_file: output/test11_4.out
458: test:
459: suffix: 17
460: nsize: 2
461: args: -resid -bv_orthog_block svqb -bv_type {{vecs contiguous svec mat}shared output}
462: output_file: output/test11_6.out
464: test:
465: suffix: 17_cuda
466: nsize: 2
467: args: -resid -bv_orthog_block svqb -bv_type svec -vec_type cuda
468: requires: cuda
469: output_file: output/test11_6.out
471: test:
472: suffix: 18
473: nsize: 2
474: args: -resid -withb -bv_orthog_block svqb -bv_type {{vecs contiguous svec mat}shared output}
475: output_file: output/test11_9.out
477: test:
478: suffix: 18_cuda
479: nsize: 2
480: args: -resid -withb -bv_orthog_block svqb -bv_type svec -vec_type cuda
481: requires: cuda
482: output_file: output/test11_9.out
484: TEST*/