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