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