Actual source code: ex118.c
petsc-3.4.5 2014-06-29
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)
14: PetscReal *work,tols[2];
15: PetscInt i,j;
16: PetscBLASInt n,il=1,iu=5,*iblock,*isplit,*iwork,nevs,*ifail,cklvl=2;
17: PetscMPIInt size;
18: PetscBool flg;
19: Vec *evecs;
20: PetscScalar *evecs_array,*D,*E,*evals;
21: Mat T;
22: #if !defined(PETSC_MISSING_LAPACK_DSTEBZ)
23: PetscReal vl=0.0,vu=4.0,tol=1.e-10;
24: PetscBLASInt nsplit,info;
25: #endif
26: #endif
28: PetscInitialize(&argc,&args,(char*)0,help);
29: #if defined(PETSC_USE_COMPLEX)
30: SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
31: #else
32: MPI_Comm_size(PETSC_COMM_WORLD,&size);
33: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
35: n = 100;
36: nevs = iu - il;
37: PetscMalloc((3*n+1)*sizeof(PetscScalar),&D);
38: E = D + n;
39: evals = E + n;
40: PetscMalloc((5*n+1)*sizeof(PetscReal),&work);
41: PetscMalloc((3*n+1)*sizeof(PetscBLASInt),&iwork);
42: PetscMalloc((3*n+1)*sizeof(PetscBLASInt),&iblock);
43: isplit = iblock + n;
45: /* Set symmetric tridiagonal matrix */
46: for (i=0; i<n; i++) {
47: D[i] = 2.0;
48: E[i] = 1.0;
49: }
51: /* Solve eigenvalue problem: A*evec = eval*evec */
52: #if defined(PETSC_MISSING_LAPACK_STEBZ)
53: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
54: #else
55: printf(" LAPACKstebz_: compute %d eigenvalues...\n",nevs);
56: LAPACKstebz_("I","E",&n,&vl,&vu,&il,&iu,&tol,(PetscReal*)D,(PetscReal*)E,&nevs,&nsplit,(PetscReal*)evals,iblock,isplit,work,iwork,&info);
57: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstebz_ fails. info %d",info);
58: #endif
60: printf(" LAPACKstein_: compute %d found eigenvectors...\n",nevs);
61: PetscMalloc(n*nevs*sizeof(PetscScalar),&evecs_array);
62: PetscMalloc(nevs*sizeof(PetscBLASInt),&ifail);
63: #if defined(PETSC_MISSING_LAPACK_STEIN)
64: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
65: #else
66: LAPACKstein_(&n,(PetscReal*)D,(PetscReal*)E,&nevs,(PetscReal*)evals,iblock,isplit,evecs_array,&n,work,iwork,ifail,&info);
67: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstein_ fails. info %d",info);
68: #endif
69: /* View evals */
70: PetscOptionsHasName(NULL, "-eig_view", &flg);
71: if (flg) {
72: PetscPrintf(PETSC_COMM_SELF," %d evals: \n",nevs);
73: for (i=0; i<nevs; i++) PetscPrintf(PETSC_COMM_SELF,"%d %G\n",i,evals[i]);
74: }
76: /* Check residuals and orthogonality */
77: MatCreate(PETSC_COMM_SELF,&T);
78: MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n);
79: MatSetType(T,MATSBAIJ);
80: MatSetFromOptions(T);
81: MatSetUp(T);
82: for (i=0; i<n; i++) {
83: MatSetValues(T,1,&i,1,&i,&D[i],INSERT_VALUES);
84: if (i != n-1) {
85: j = i+1;
86: MatSetValues(T,1,&i,1,&j,&E[i],INSERT_VALUES);
87: }
88: }
89: MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);
92: PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
93: for (i=0; i<nevs; i++) {
94: VecCreate(PETSC_COMM_SELF,&evecs[i]);
95: VecSetSizes(evecs[i],PETSC_DECIDE,n);
96: VecSetFromOptions(evecs[i]);
97: VecPlaceArray(evecs[i],evecs_array+i*n);
98: }
100: tols[0] = 1.e-8; tols[1] = 1.e-8;
101: CkEigenSolutions(cklvl,T,il-1,iu-1,evals,evecs,tols);
103: for (i=0; i<nevs; i++) {
104: VecResetArray(evecs[i]);
105: }
107: /* free space */
109: MatDestroy(&T);
111: for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
112: PetscFree(evecs);
113: PetscFree(D);
114: PetscFree(work);
115: PetscFree(iwork);
116: PetscFree(iblock);
117: PetscFree(evecs_array);
118: PetscFree(ifail);
119: PetscFinalize();
120: #endif
121: return 0;
122: }
123: /*------------------------------------------------
124: Check the accuracy of the eigen solution
125: ----------------------------------------------- */
126: /*
127: input:
128: cklvl - check level:
129: 1: check residual
130: 2: 1 and check B-orthogonality locally
131: A - matrix
132: il,iu - lower and upper index bound of eigenvalues
133: eval, evec - eigenvalues and eigenvectors stored in this process
134: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
135: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
136: */
137: #undef DEBUG_CkEigenSolutions
140: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscScalar *eval,Vec *evec,PetscReal *tols)
141: {
142: PetscInt ierr,i,j,nev;
143: Vec vt1,vt2; /* tmp vectors */
144: PetscReal norm,norm_max;
145: PetscScalar dot,tmp;
146: PetscReal dot_max;
149: nev = iu - il;
150: if (nev <= 0) return(0);
152: VecDuplicate(evec[0],&vt1);
153: VecDuplicate(evec[0],&vt2);
155: switch (cklvl) {
156: case 2:
157: dot_max = 0.0;
158: for (i = il; i<iu; i++) {
159: VecCopy(evec[i], vt1);
160: for (j=il; j<iu; j++) {
161: VecDot(evec[j],vt1,&dot);
162: if (j == i) {
163: dot = PetscAbsScalar(dot - 1.0);
164: } else {
165: dot = PetscAbsScalar(dot);
166: }
167: if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot);
168: #if defined(DEBUG_CkEigenSolutions)
169: if (dot > tols[1]) {
170: VecNorm(evec[i],NORM_INFINITY,&norm);
171: PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
172: }
173: #endif
174: }
175: }
176: PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %G\n",dot_max);
178: case 1:
179: norm_max = 0.0;
180: for (i = il; i< iu; i++) {
181: MatMult(A, evec[i], vt1);
182: VecCopy(evec[i], vt2);
183: tmp = -eval[i];
184: VecAXPY(vt1,tmp,vt2);
185: VecNorm(vt1, NORM_INFINITY, &norm);
186: norm = PetscAbsScalar(norm);
187: if (norm > norm_max) norm_max = norm;
188: #if defined(DEBUG_CkEigenSolutions)
189: /* sniff, and bark if necessary */
190: if (norm > tols[0]) {
191: printf(" residual violation: %d, resi: %g\n",i, norm);
192: }
193: #endif
194: }
195: PetscPrintf(PETSC_COMM_SELF," max_resi: %G\n", norm_max);
196: break;
197: default:
198: PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
199: }
201: VecDestroy(&vt2);
202: VecDestroy(&vt1);
203: return(0);
204: }