Actual source code: ex118.c
petsc-3.6.4 2016-04-12
1: static char help[] = "Test LAPACK routine DSTEBZ() and DTEIN(). \n\n";
3: #include <petscmat.h>
4: #include <petscblaslapack.h>
6: extern PetscErrorCode CkEigenSolutions(PetscInt,Mat,PetscInt,PetscInt,PetscScalar*,Vec*,PetscReal*);
10: int main(int argc,char **args)
11: {
12: #if defined(PETSC_USE_COMPLEX) || defined(PETSC_MISSING_LAPACK_DSTEBZ) || defined(PETSC_MISSING_LAPACK_STEIN)
13: PetscInitialize(&argc,&args,(char*)0,help);
14: SETERRQ(PETSC_COMM_WORLD,1,"This example requires LAPACK routines dstebz and stien and real numbers");
15: #else
17: PetscReal *work,tols[2];
18: PetscInt i,j;
19: PetscBLASInt n,il=1,iu=5,*iblock,*isplit,*iwork,nevs,*ifail,cklvl=2;
20: PetscMPIInt size;
21: PetscBool flg;
22: Vec *evecs;
23: PetscScalar *evecs_array,*D,*E,*evals;
24: Mat T;
25: PetscReal vl=0.0,vu=4.0,tol=1.e-10;
26: PetscBLASInt nsplit,info;
29: PetscInitialize(&argc,&args,(char*)0,help);
30: MPI_Comm_size(PETSC_COMM_WORLD,&size);
31: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
33: n = 100;
34: nevs = iu - il;
35: PetscMalloc1(3*n+1,&D);
36: E = D + n;
37: evals = E + n;
38: PetscMalloc1(5*n+1,&work);
39: PetscMalloc1(3*n+1,&iwork);
40: PetscMalloc1(3*n+1,&iblock);
41: isplit = iblock + n;
43: /* Set symmetric tridiagonal matrix */
44: for (i=0; i<n; i++) {
45: D[i] = 2.0;
46: E[i] = 1.0;
47: }
49: /* Solve eigenvalue problem: A*evec = eval*evec */
50: PetscPrintf(PETSC_COMM_SELF," LAPACKstebz_: compute %d eigenvalues...\n",nevs);
51: LAPACKstebz_("I","E",&n,&vl,&vu,&il,&iu,&tol,(PetscReal*)D,(PetscReal*)E,&nevs,&nsplit,(PetscReal*)evals,iblock,isplit,work,iwork,&info);
52: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstebz_ fails. info %d",info);
54: PetscPrintf(PETSC_COMM_SELF," LAPACKstein_: compute %d found eigenvectors...\n",nevs);
55: PetscMalloc1(n*nevs,&evecs_array);
56: PetscMalloc1(nevs,&ifail);
57: LAPACKstein_(&n,(PetscReal*)D,(PetscReal*)E,&nevs,(PetscReal*)evals,iblock,isplit,evecs_array,&n,work,iwork,ifail,&info);
58: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstein_ fails. info %d",info);
59: /* View evals */
60: PetscOptionsHasName(NULL, "-eig_view", &flg);
61: if (flg) {
62: PetscPrintf(PETSC_COMM_SELF," %d evals: \n",nevs);
63: for (i=0; i<nevs; i++) {PetscPrintf(PETSC_COMM_SELF,"%D %g\n",i,(double)evals[i]);}
64: }
66: /* Check residuals and orthogonality */
67: MatCreate(PETSC_COMM_SELF,&T);
68: MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n);
69: MatSetType(T,MATSBAIJ);
70: MatSetFromOptions(T);
71: MatSetUp(T);
72: for (i=0; i<n; i++) {
73: MatSetValues(T,1,&i,1,&i,&D[i],INSERT_VALUES);
74: if (i != n-1) {
75: j = i+1;
76: MatSetValues(T,1,&i,1,&j,&E[i],INSERT_VALUES);
77: }
78: }
79: MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);
82: PetscMalloc1(nevs+1,&evecs);
83: for (i=0; i<nevs; i++) {
84: VecCreate(PETSC_COMM_SELF,&evecs[i]);
85: VecSetSizes(evecs[i],PETSC_DECIDE,n);
86: VecSetFromOptions(evecs[i]);
87: VecPlaceArray(evecs[i],evecs_array+i*n);
88: }
90: tols[0] = 1.e-8; tols[1] = 1.e-8;
91: CkEigenSolutions(cklvl,T,il-1,iu-1,evals,evecs,tols);
93: for (i=0; i<nevs; i++) {
94: VecResetArray(evecs[i]);
95: }
97: /* free space */
99: MatDestroy(&T);
101: for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
102: PetscFree(evecs);
103: PetscFree(D);
104: PetscFree(work);
105: PetscFree(iwork);
106: PetscFree(iblock);
107: PetscFree(evecs_array);
108: PetscFree(ifail);
109: PetscFinalize();
110: return 0;
111: #endif
112: }
113: /*------------------------------------------------
114: Check the accuracy of the eigen solution
115: ----------------------------------------------- */
116: /*
117: input:
118: cklvl - check level:
119: 1: check residual
120: 2: 1 and check B-orthogonality locally
121: A - matrix
122: il,iu - lower and upper index bound of eigenvalues
123: eval, evec - eigenvalues and eigenvectors stored in this process
124: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
125: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
126: */
127: #undef DEBUG_CkEigenSolutions
130: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscScalar *eval,Vec *evec,PetscReal *tols)
131: {
132: PetscInt ierr,i,j,nev;
133: Vec vt1,vt2; /* tmp vectors */
134: PetscReal norm,norm_max;
135: PetscScalar dot,tmp;
136: PetscReal dot_max;
139: nev = iu - il;
140: if (nev <= 0) return(0);
142: VecDuplicate(evec[0],&vt1);
143: VecDuplicate(evec[0],&vt2);
145: switch (cklvl) {
146: case 2:
147: dot_max = 0.0;
148: for (i = il; i<iu; i++) {
149: VecCopy(evec[i], vt1);
150: for (j=il; j<iu; j++) {
151: VecDot(evec[j],vt1,&dot);
152: if (j == i) {
153: dot = PetscAbsScalar(dot - 1.0);
154: } else {
155: dot = PetscAbsScalar(dot);
156: }
157: if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot);
158: #if defined(DEBUG_CkEigenSolutions)
159: if (dot > tols[1]) {
160: VecNorm(evec[i],NORM_INFINITY,&norm);
161: PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)dot,(double)norm);
162: }
163: #endif
164: }
165: }
166: PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);
168: case 1:
169: norm_max = 0.0;
170: for (i = il; i< iu; i++) {
171: MatMult(A, evec[i], vt1);
172: VecCopy(evec[i], vt2);
173: tmp = -eval[i];
174: VecAXPY(vt1,tmp,vt2);
175: VecNorm(vt1, NORM_INFINITY, &norm);
176: norm = PetscAbsReal(norm);
177: if (norm > norm_max) norm_max = norm;
178: #if defined(DEBUG_CkEigenSolutions)
179: if (norm > tols[0]) {
180: PetscPrintf(PETSC_COMM_SELF," residual violation: %d, resi: %g\n",i, norm);
181: }
182: #endif
183: }
184: PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max);
185: break;
186: default:
187: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
188: }
190: VecDestroy(&vt2);
191: VecDestroy(&vt1);
192: return(0);
193: }