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)
12: PetscInitialize(&argc, &args, (char *)0, help);
13: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, "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;
27: PetscInitialize(&argc, &args, (char *)0, help);
28: MPI_Comm_size(PETSC_COMM_WORLD, &size);
31: n = 100;
32: nevs = iu - il;
33: PetscMalloc1(3 * n + 1, &D);
34: E = D + n;
35: evals = E + n;
36: PetscMalloc1(5 * n + 1, &work);
37: PetscMalloc1(3 * n + 1, &iwork);
38: PetscMalloc1(3 * n + 1, &iblock);
39: isplit = iblock + n;
41: /* Set symmetric tridiagonal matrix */
42: for (i = 0; i < n; i++) {
43: D[i] = 2.0;
44: E[i] = 1.0;
45: }
47: /* Solve eigenvalue problem: A*evec = eval*evec */
48: PetscPrintf(PETSC_COMM_SELF, " LAPACKstebz_: compute %d eigenvalues...\n", nevs);
49: LAPACKstebz_("I", "E", &n, &vl, &vu, &il, &iu, &tol, (PetscReal *)D, (PetscReal *)E, &nevs, &nsplit, (PetscReal *)evals, iblock, isplit, work, iwork, &info);
52: PetscPrintf(PETSC_COMM_SELF, " LAPACKstein_: compute %d found eigenvectors...\n", nevs);
53: PetscMalloc1(n * nevs, &evecs_array);
54: PetscMalloc1(nevs, &ifail);
55: LAPACKstein_(&n, (PetscReal *)D, (PetscReal *)E, &nevs, (PetscReal *)evals, iblock, isplit, evecs_array, &n, work, iwork, ifail, &info);
57: /* View evals */
58: PetscOptionsHasName(NULL, NULL, "-eig_view", &flg);
59: if (flg) {
60: PetscPrintf(PETSC_COMM_SELF, " %d evals: \n", nevs);
61: for (i = 0; i < nevs; i++) PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %g\n", i, (double)evals[i]);
62: }
64: /* Check residuals and orthogonality */
65: MatCreate(PETSC_COMM_SELF, &T);
66: MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, n, n);
67: MatSetType(T, MATSBAIJ);
68: MatSetFromOptions(T);
69: MatSetUp(T);
70: for (i = 0; i < n; i++) {
71: MatSetValues(T, 1, &i, 1, &i, &D[i], INSERT_VALUES);
72: if (i != n - 1) {
73: j = i + 1;
74: MatSetValues(T, 1, &i, 1, &j, &E[i], INSERT_VALUES);
75: }
76: }
77: MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY);
80: PetscMalloc1(nevs + 1, &evecs);
81: for (i = 0; i < nevs; i++) {
82: VecCreate(PETSC_COMM_SELF, &evecs[i]);
83: VecSetSizes(evecs[i], PETSC_DECIDE, n);
84: VecSetFromOptions(evecs[i]);
85: VecPlaceArray(evecs[i], evecs_array + i * n);
86: }
88: tols[0] = 1.e-8;
89: tols[1] = 1.e-8;
90: CkEigenSolutions(cklvl, T, il - 1, iu - 1, evals, evecs, tols);
92: for (i = 0; i < nevs; i++) VecResetArray(evecs[i]);
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 0;
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;
133: nev = iu - il;
134: if (nev <= 0) return 0;
136: VecDuplicate(evec[0], &vt1);
137: VecDuplicate(evec[0], &vt2);
139: switch (cklvl) {
140: case 2:
141: dot_max = 0.0;
142: for (i = il; i < iu; i++) {
143: VecCopy(evec[i], vt1);
144: for (j = il; j < iu; j++) {
145: VecDot(evec[j], vt1, &dot);
146: if (j == i) {
147: dot = PetscAbsScalar(dot - (PetscScalar)1.0);
148: } else {
149: dot = PetscAbsScalar(dot);
150: }
151: if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot);
152: #if defined(DEBUG_CkEigenSolutions)
153: if (dot > tols[1]) {
154: VecNorm(evec[i], NORM_INFINITY, &norm);
155: PetscPrintf(PETSC_COMM_SELF, "|delta(%d,%d)|: %g, norm: %d\n", i, j, (double)dot, (double)norm);
156: }
157: #endif
158: }
159: }
160: PetscPrintf(PETSC_COMM_SELF, " max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max);
162: case 1:
163: norm_max = 0.0;
164: for (i = il; i < iu; i++) {
165: MatMult(A, evec[i], vt1);
166: VecCopy(evec[i], vt2);
167: tmp = -eval[i];
168: VecAXPY(vt1, tmp, vt2);
169: VecNorm(vt1, NORM_INFINITY, &norm);
170: norm = PetscAbsReal(norm);
171: if (norm > norm_max) norm_max = norm;
172: #if defined(DEBUG_CkEigenSolutions)
173: if (norm > tols[0]) PetscPrintf(PETSC_COMM_SELF, " residual violation: %d, resi: %g\n", i, norm);
174: #endif
175: }
176: PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max);
177: break;
178: default:
179: PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%d is not supported \n", cklvl);
180: }
182: VecDestroy(&vt2);
183: VecDestroy(&vt1);
184: return 0;
185: }