1 static char help[] = "Test LAPACK routine DSTEBZ() and DTEIN(). \n\n";
2
3 #include <petscmat.h>
4 #include <petscblaslapack.h>
5
6 extern PetscErrorCode CkEigenSolutions(PetscInt, Mat, PetscInt, PetscInt, PetscScalar *, Vec *, PetscReal *);
7
main(int argc,char ** args)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;
25
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!");
30
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;
40
41 /* Set symmetric tridiagonal matrix */
42 for (i = 0; i < n; i++) {
43 D[i] = 2.0;
44 E[i] = 1.0;
45 }
46
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);
51
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 }
63
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));
79
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 }
87
88 tols[0] = 1.e-8;
89 tols[1] = 1.e-8;
90 PetscCall(CkEigenSolutions(cklvl, T, il - 1, iu - 1, evals, evecs, tols));
91
92 for (i = 0; i < nevs; i++) PetscCall(VecResetArray(evecs[i]));
93
94 /* free space */
95
96 PetscCall(MatDestroy(&T));
97
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
CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscScalar * eval,Vec * evec,PetscReal * tols)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;
132
133 PetscFunctionBegin;
134 nev = iu - il;
135 if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);
136
137 PetscCall(VecDuplicate(evec[0], &vt1));
138 PetscCall(VecDuplicate(evec[0], &vt2));
139
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));
162
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 }
182
183 PetscCall(VecDestroy(&vt2));
184 PetscCall(VecDestroy(&vt1));
185 PetscFunctionReturn(PETSC_SUCCESS);
186 }
187