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