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: PetscFunctionBeginUser;
12: PetscCall(PetscInitialize(&argc, &args, NULL, 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;
26: PetscFunctionBeginUser;
27: PetscCall(PetscInitialize(&argc, &args, NULL, help));
28: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
29: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
31: n = 100;
32: nevs = iu - il;
33: PetscCall(PetscMalloc1(3 * n + 1, &D));
34: E = D + n;
35: evals = E + n;
36: PetscCall(PetscMalloc1(5 * n + 1, &work));
37: PetscCall(PetscMalloc1(3 * n + 1, &iwork));
38: PetscCall(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: PetscCall(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);
50: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_USER, "LAPACKstebz_ fails. info %d", info);
52: PetscCall(PetscPrintf(PETSC_COMM_SELF, " LAPACKstein_: compute %d found eigenvectors...\n", nevs));
53: PetscCall(PetscMalloc1(n * nevs, &evecs_array));
54: PetscCall(PetscMalloc1(nevs, &ifail));
55: LAPACKstein_(&n, (PetscReal *)D, (PetscReal *)E, &nevs, (PetscReal *)evals, iblock, isplit, evecs_array, &n, work, iwork, ifail, &info);
56: PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_USER, "LAPACKstein_ fails. info %d", info);
57: /* View evals */
58: PetscCall(PetscOptionsHasName(NULL, NULL, "-eig_view", &flg));
59: if (flg) {
60: PetscCall(PetscPrintf(PETSC_COMM_SELF, " %d evals: \n", nevs));
61: for (i = 0; i < nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %g\n", i, (double)evals[i]));
62: }
64: /* Check residuals and orthogonality */
65: PetscCall(MatCreate(PETSC_COMM_SELF, &T));
66: PetscCall(MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, n, n));
67: PetscCall(MatSetType(T, MATSBAIJ));
68: PetscCall(MatSetFromOptions(T));
69: PetscCall(MatSetUp(T));
70: for (i = 0; i < n; i++) {
71: PetscCall(MatSetValues(T, 1, &i, 1, &i, &D[i], INSERT_VALUES));
72: if (i != n - 1) {
73: j = i + 1;
74: PetscCall(MatSetValues(T, 1, &i, 1, &j, &E[i], INSERT_VALUES));
75: }
76: }
77: PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
78: PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));
80: PetscCall(PetscMalloc1(nevs + 1, &evecs));
81: for (i = 0; i < nevs; i++) {
82: PetscCall(VecCreate(PETSC_COMM_SELF, &evecs[i]));
83: PetscCall(VecSetSizes(evecs[i], PETSC_DECIDE, n));
84: PetscCall(VecSetFromOptions(evecs[i]));
85: PetscCall(VecPlaceArray(evecs[i], evecs_array + i * n));
86: }
88: tols[0] = 1.e-8;
89: tols[1] = 1.e-8;
90: PetscCall(CkEigenSolutions(cklvl, T, il - 1, iu - 1, evals, evecs, tols));
92: for (i = 0; i < nevs; i++) PetscCall(VecResetArray(evecs[i]));
94: /* free space */
96: PetscCall(MatDestroy(&T));
98: for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i]));
99: PetscCall(PetscFree(evecs));
100: PetscCall(PetscFree(D));
101: PetscCall(PetscFree(work));
102: PetscCall(PetscFree(iwork));
103: PetscCall(PetscFree(iblock));
104: PetscCall(PetscFree(evecs_array));
105: PetscCall(PetscFree(ifail));
106: PetscCall(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: PetscFunctionBegin;
134: nev = iu - il;
135: if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);
137: PetscCall(VecDuplicate(evec[0], &vt1));
138: PetscCall(VecDuplicate(evec[0], &vt2));
140: switch (cklvl) {
141: case 2:
142: dot_max = 0.0;
143: for (i = il; i < iu; i++) {
144: PetscCall(VecCopy(evec[i], vt1));
145: for (j = il; j < iu; j++) {
146: PetscCall(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: PetscCall(VecNorm(evec[i], NORM_INFINITY, &norm));
156: PetscCall(PetscPrintf(PETSC_COMM_SELF, "|delta(%d,%d)|: %g, norm: %d\n", i, j, (double)dot, (double)norm));
157: }
158: #endif
159: }
160: }
161: PetscCall(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: PetscCall(MatMult(A, evec[i], vt1));
167: PetscCall(VecCopy(evec[i], vt2));
168: tmp = -eval[i];
169: PetscCall(VecAXPY(vt1, tmp, vt2));
170: PetscCall(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]) PetscCall(PetscPrintf(PETSC_COMM_SELF, " residual violation: %d, resi: %g\n", i, norm));
175: #endif
176: }
177: PetscCall(PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max));
178: break;
179: default:
180: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%d is not supported \n", cklvl));
181: }
183: PetscCall(VecDestroy(&vt2));
184: PetscCall(VecDestroy(&vt1));
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }