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 PetscDefined(USE_COMPLEX) || PetscDefined(MISSING_LAPACK_STEBZ) || PetscDefined(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;
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: PetscCallLAPACKInfo("LAPACKstebz", LAPACKstebz_("I", "E", &n, &vl, &vu, &il, &iu, &tol, (PetscReal *)D, (PetscReal *)E, &nevs, &nsplit, (PetscReal *)evals, iblock, isplit, work, iwork, &info));
50: PetscCall(PetscPrintf(PETSC_COMM_SELF, " LAPACKstein_: compute %d found eigenvectors...\n", nevs));
51: PetscCall(PetscMalloc1(n * nevs, &evecs_array));
52: PetscCall(PetscMalloc1(nevs, &ifail));
53: PetscCallLAPACKInfo("LAPACKstein", LAPACKstein_(&n, (PetscReal *)D, (PetscReal *)E, &nevs, (PetscReal *)evals, iblock, isplit, evecs_array, &n, work, iwork, ifail, &info));
55: /* View evals */
56: PetscCall(PetscOptionsHasName(NULL, NULL, "-eig_view", &flg));
57: if (flg) {
58: PetscCall(PetscPrintf(PETSC_COMM_SELF, " %d evals: \n", nevs));
59: for (i = 0; i < nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %g\n", i, (double)evals[i]));
60: }
62: /* Check residuals and orthogonality */
63: PetscCall(MatCreate(PETSC_COMM_SELF, &T));
64: PetscCall(MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, n, n));
65: PetscCall(MatSetType(T, MATSBAIJ));
66: PetscCall(MatSetFromOptions(T));
67: PetscCall(MatSetUp(T));
68: for (i = 0; i < n; i++) {
69: PetscCall(MatSetValues(T, 1, &i, 1, &i, &D[i], INSERT_VALUES));
70: if (i != n - 1) {
71: j = i + 1;
72: PetscCall(MatSetValues(T, 1, &i, 1, &j, &E[i], INSERT_VALUES));
73: }
74: }
75: PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY));
76: PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY));
78: PetscCall(PetscMalloc1(nevs + 1, &evecs));
79: for (i = 0; i < nevs; i++) {
80: PetscCall(VecCreate(PETSC_COMM_SELF, &evecs[i]));
81: PetscCall(VecSetSizes(evecs[i], PETSC_DECIDE, n));
82: PetscCall(VecSetFromOptions(evecs[i]));
83: PetscCall(VecPlaceArray(evecs[i], evecs_array + i * n));
84: }
86: tols[0] = 1.e-8;
87: tols[1] = 1.e-8;
88: PetscCall(CkEigenSolutions(cklvl, T, il - 1, iu - 1, evals, evecs, tols));
90: for (i = 0; i < nevs; i++) PetscCall(VecResetArray(evecs[i]));
92: /* free space */
94: PetscCall(MatDestroy(&T));
96: for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i]));
97: PetscCall(PetscFree(evecs));
98: PetscCall(PetscFree(D));
99: PetscCall(PetscFree(work));
100: PetscCall(PetscFree(iwork));
101: PetscCall(PetscFree(iblock));
102: PetscCall(PetscFree(evecs_array));
103: PetscCall(PetscFree(ifail));
104: PetscCall(PetscFinalize());
105: return 0;
106: #endif
107: }
108: /*------------------------------------------------
109: Check the accuracy of the eigen solution
110: ----------------------------------------------- */
111: /*
112: input:
113: cklvl - check level:
114: 1: check residual
115: 2: 1 and check B-orthogonality locally
116: A - matrix
117: il,iu - lower and upper index bound of eigenvalues
118: eval, evec - eigenvalues and eigenvectors stored in this process
119: tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
120: tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
121: */
122: #undef DEBUG_CkEigenSolutions
123: PetscErrorCode CkEigenSolutions(PetscInt cklvl, Mat A, PetscInt il, PetscInt iu, PetscScalar *eval, Vec *evec, PetscReal *tols)
124: {
125: PetscInt i, j, nev;
126: Vec vt1, vt2; /* tmp vectors */
127: PetscReal norm, norm_max;
128: PetscScalar dot, tmp;
129: PetscReal dot_max;
131: PetscFunctionBegin;
132: nev = iu - il;
133: if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);
135: PetscCall(VecDuplicate(evec[0], &vt1));
136: PetscCall(VecDuplicate(evec[0], &vt2));
138: switch (cklvl) {
139: case 2:
140: dot_max = 0.0;
141: for (i = il; i < iu; i++) {
142: PetscCall(VecCopy(evec[i], vt1));
143: for (j = il; j < iu; j++) {
144: PetscCall(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: PetscCall(VecNorm(evec[i], NORM_INFINITY, &norm));
154: PetscCall(PetscPrintf(PETSC_COMM_SELF, "|delta(%d,%d)|: %g, norm: %d\n", i, j, (double)dot, (double)norm));
155: }
156: #endif
157: }
158: }
159: PetscCall(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: PetscCall(MatMult(A, evec[i], vt1));
165: PetscCall(VecCopy(evec[i], vt2));
166: tmp = -eval[i];
167: PetscCall(VecAXPY(vt1, tmp, vt2));
168: PetscCall(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]) PetscCall(PetscPrintf(PETSC_COMM_SELF, " residual violation: %d, resi: %g\n", i, norm));
173: #endif
174: }
175: PetscCall(PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max));
176: break;
177: default:
178: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%d is not supported \n", cklvl));
179: }
181: PetscCall(VecDestroy(&vt2));
182: PetscCall(VecDestroy(&vt1));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }