1 static char help[] = "Test LAPACK routine DSYEV() or DSYEVX(). \n\
2 Reads PETSc matrix A \n\
3 then computes selected eigenvalues, and optionally, eigenvectors of \n\
4 a real generalized symmetric-definite eigenproblem \n\
5 A*x = lambda*x \n\
6 Input parameters include\n\
7 -f <input_file> : file to load\n\
8 e.g. ./ex116 -f $DATAFILESPATH/matrices/small \n\n";
9
10 #include <petscmat.h>
11 #include <petscblaslapack.h>
12
13 extern PetscErrorCode CkEigenSolutions(PetscInt, Mat, PetscInt, PetscInt, PetscReal *, Vec *, PetscReal *);
14
main(int argc,char ** args)15 int main(int argc, char **args)
16 {
17 Mat A, A_dense;
18 Vec *evecs;
19 PetscViewer fd; /* viewer */
20 char file[1][PETSC_MAX_PATH_LEN]; /* input file name */
21 PetscBool flg, TestSYEVX = PETSC_TRUE;
22 PetscBool isSymmetric;
23 PetscScalar *arrayA, *evecs_array, *work, *evals;
24 PetscMPIInt size;
25 PetscInt m, n, i, cklvl = 2;
26 PetscBLASInt nevs, il, iu, in;
27 PetscReal vl, vu, abstol = 1.e-8;
28 PetscBLASInt *iwork, *ifail, lwork, lierr, bn;
29 PetscReal tols[2];
30
31 PetscFunctionBeginUser;
32 PetscCall(PetscInitialize(&argc, &args, NULL, help));
33 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
34 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
35
36 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_syev", &flg));
37 if (flg) TestSYEVX = PETSC_FALSE;
38
39 /* Determine files from which we read the two matrices */
40 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg));
41
42 /* Load matrix A */
43 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &fd));
44 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
45 PetscCall(MatSetType(A, MATSEQAIJ));
46 PetscCall(MatLoad(A, fd));
47 PetscCall(PetscViewerDestroy(&fd));
48 PetscCall(MatGetSize(A, &m, &n));
49
50 /* Check whether A is symmetric */
51 PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetry", &flg));
52 if (flg) {
53 Mat Trans;
54 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Trans));
55 PetscCall(MatEqual(A, Trans, &isSymmetric));
56 PetscCheck(isSymmetric, PETSC_COMM_SELF, PETSC_ERR_USER, "A must be symmetric");
57 PetscCall(MatDestroy(&Trans));
58 }
59
60 /* Solve eigenvalue problem: A_dense*x = lambda*B*x */
61 /*==================================================*/
62 /* Convert aij matrix to MatSeqDense for LAPACK */
63 PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense));
64
65 PetscCall(PetscBLASIntCast(8 * n, &lwork));
66 PetscCall(PetscBLASIntCast(n, &bn));
67 PetscCall(PetscMalloc1(n, &evals));
68 PetscCall(PetscMalloc1(lwork, &work));
69 PetscCall(MatDenseGetArray(A_dense, &arrayA));
70
71 if (!TestSYEVX) { /* test syev() */
72 PetscCall(PetscPrintf(PETSC_COMM_SELF, " LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n", m));
73 LAPACKsyev_("V", "U", &bn, arrayA, &bn, evals, work, &lwork, &lierr);
74 evecs_array = arrayA;
75 PetscCall(PetscBLASIntCast(m, &nevs));
76 il = 1;
77 PetscCall(PetscBLASIntCast(m, &iu));
78 } else { /* test syevx() */
79 il = 1;
80 PetscCall(PetscBLASIntCast(0.2 * m, &iu));
81 PetscCall(PetscBLASIntCast(n, &in));
82 PetscCall(PetscPrintf(PETSC_COMM_SELF, " LAPACKsyevx: compute %" PetscBLASInt_FMT " to %" PetscBLASInt_FMT "-th eigensolutions...\n", il, iu));
83 PetscCall(PetscMalloc1(m * n + 1, &evecs_array));
84 PetscCall(PetscMalloc1(6 * n + 1, &iwork));
85 ifail = iwork + 5 * n;
86
87 /* in the case "I", vl and vu are not referenced */
88 vl = 0.0;
89 vu = 8.0;
90 LAPACKsyevx_("V", "I", "U", &bn, arrayA, &bn, &vl, &vu, &il, &iu, &abstol, &nevs, evals, evecs_array, &in, work, &lwork, iwork, ifail, &lierr);
91 PetscCall(PetscFree(iwork));
92 }
93 PetscCall(MatDenseRestoreArray(A_dense, &arrayA));
94 PetscCheck(nevs > 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "nev=%" PetscBLASInt_FMT ", no eigensolution has found", nevs);
95
96 /* View eigenvalues */
97 PetscCall(PetscOptionsHasName(NULL, NULL, "-eig_view", &flg));
98 if (flg) {
99 PetscCall(PetscPrintf(PETSC_COMM_SELF, " %" PetscBLASInt_FMT " evals: \n", nevs));
100 for (i = 0; i < nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %g\n", (PetscInt)(i + il), (double)evals[i]));
101 }
102
103 /* Check residuals and orthogonality */
104 PetscCall(PetscMalloc1(nevs + 1, &evecs));
105 for (i = 0; i < nevs; i++) {
106 PetscCall(VecCreate(PETSC_COMM_SELF, &evecs[i]));
107 PetscCall(VecSetSizes(evecs[i], PETSC_DECIDE, n));
108 PetscCall(VecSetFromOptions(evecs[i]));
109 PetscCall(VecPlaceArray(evecs[i], evecs_array + i * n));
110 }
111
112 tols[0] = tols[1] = PETSC_SQRT_MACHINE_EPSILON;
113 PetscCall(CkEigenSolutions(cklvl, A, il - 1, iu - 1, evals, evecs, tols));
114
115 /* Free work space. */
116 for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i]));
117 PetscCall(PetscFree(evecs));
118 PetscCall(MatDestroy(&A_dense));
119 PetscCall(PetscFree(work));
120 if (TestSYEVX) PetscCall(PetscFree(evecs_array));
121
122 /* Compute SVD: A_dense = U*SIGMA*transpose(V),
123 JOBU=JOBV='S': the first min(m,n) columns of U and V are returned in the arrayU and arrayV; */
124 /*==============================================================================================*/
125 {
126 /* Convert aij matrix to MatSeqDense for LAPACK */
127 PetscScalar *arrayU, *arrayVT, *arrayErr, alpha = 1.0, beta = -1.0;
128 Mat Err;
129 PetscBLASInt minMN, maxMN, im, in;
130 PetscInt j;
131 PetscReal norm;
132
133 PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense));
134
135 PetscCall(PetscBLASIntCast(PetscMin(m, n), &minMN));
136 PetscCall(PetscBLASIntCast(PetscMax(m, n), &maxMN));
137 PetscCall(PetscBLASIntCast(5 * minMN + maxMN, &lwork));
138 PetscCall(PetscMalloc4(m * minMN, &arrayU, m * minMN, &arrayVT, m * minMN, &arrayErr, lwork, &work));
139
140 /* Create matrix Err for checking error */
141 PetscCall(MatCreate(PETSC_COMM_WORLD, &Err));
142 PetscCall(MatSetSizes(Err, PETSC_DECIDE, PETSC_DECIDE, m, minMN));
143 PetscCall(MatSetType(Err, MATSEQDENSE));
144 PetscCall(MatSeqDenseSetPreallocation(Err, (PetscScalar *)arrayErr));
145
146 /* Save A to arrayErr for checking accuracy later. arrayA will be destroyed by LAPACKgesvd_() */
147 PetscCall(MatDenseGetArray(A_dense, &arrayA));
148 PetscCall(PetscArraycpy(arrayErr, arrayA, m * minMN));
149
150 PetscCall(PetscBLASIntCast(m, &im));
151 PetscCall(PetscBLASIntCast(n, &in));
152 /* Compute A = U*SIGMA*VT */
153 LAPACKgesvd_("S", "S", &im, &in, arrayA, &im, evals, arrayU, &minMN, arrayVT, &minMN, work, &lwork, &lierr);
154 PetscCall(MatDenseRestoreArray(A_dense, &arrayA));
155 if (!lierr) {
156 PetscCall(PetscPrintf(PETSC_COMM_SELF, " 1st 10 of %" PetscBLASInt_FMT " singular values: \n", minMN));
157 for (i = 0; i < 10; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %g\n", i, (double)evals[i]));
158 } else {
159 PetscCall(PetscPrintf(PETSC_COMM_SELF, "LAPACKgesvd_ fails!"));
160 }
161
162 /* Check Err = (U*Sigma*V^T - A) using BLASgemm() */
163 /* U = U*Sigma */
164 for (j = 0; j < minMN; j++) { /* U[:,j] = sigma[j]*U[:,j] */
165 for (i = 0; i < m; i++) arrayU[j * m + i] *= evals[j];
166 }
167 /* Err = U*VT - A = alpha*U*VT + beta*Err */
168 BLASgemm_("N", "N", &im, &minMN, &minMN, &alpha, arrayU, &im, arrayVT, &minMN, &beta, arrayErr, &im);
169 PetscCall(MatNorm(Err, NORM_FROBENIUS, &norm));
170 PetscCall(PetscPrintf(PETSC_COMM_SELF, " || U*Sigma*VT - A || = %g\n", (double)norm));
171
172 PetscCall(PetscFree4(arrayU, arrayVT, arrayErr, work));
173 PetscCall(PetscFree(evals));
174 PetscCall(MatDestroy(&A_dense));
175 PetscCall(MatDestroy(&Err));
176 }
177
178 PetscCall(MatDestroy(&A));
179 PetscCall(PetscFinalize());
180 return 0;
181 }
182 /*------------------------------------------------
183 Check the accuracy of the eigen solution
184 ----------------------------------------------- */
185 /*
186 input:
187 cklvl - check level:
188 1: check residual
189 2: 1 and check B-orthogonality locally
190 A - matrix
191 il,iu - lower and upper index bound of eigenvalues
192 eval, evec - eigenvalues and eigenvectors stored in this process
193 tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
194 tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
195 */
CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal * eval,Vec * evec,PetscReal * tols)196 PetscErrorCode CkEigenSolutions(PetscInt cklvl, Mat A, PetscInt il, PetscInt iu, PetscReal *eval, Vec *evec, PetscReal *tols)
197 {
198 PetscInt i, j, nev;
199 Vec vt1, vt2; /* tmp vectors */
200 PetscReal norm, tmp, dot, norm_max, dot_max;
201
202 PetscFunctionBegin;
203 nev = iu - il;
204 if (nev <= 0) PetscFunctionReturn(PETSC_SUCCESS);
205
206 /*ierr = VecView(evec[0],PETSC_VIEWER_STDOUT_WORLD);*/
207 PetscCall(VecDuplicate(evec[0], &vt1));
208 PetscCall(VecDuplicate(evec[0], &vt2));
209
210 switch (cklvl) {
211 case 2:
212 dot_max = 0.0;
213 for (i = il; i < iu; i++) {
214 PetscCall(VecCopy(evec[i], vt1));
215 for (j = il; j < iu; j++) {
216 PetscCall(VecDot(evec[j], vt1, &dot));
217 if (j == i) {
218 dot = PetscAbsScalar(dot - 1);
219 } else {
220 dot = PetscAbsScalar(dot);
221 }
222 if (dot > dot_max) dot_max = dot;
223 if (dot > tols[1]) {
224 PetscCall(VecNorm(evec[i], NORM_INFINITY, &norm));
225 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n", i, j, (double)dot, (double)norm));
226 }
227 }
228 }
229 PetscCall(PetscPrintf(PETSC_COMM_SELF, " max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max));
230
231 case 1:
232 norm_max = 0.0;
233 for (i = il; i < iu; i++) {
234 PetscCall(MatMult(A, evec[i], vt1));
235 PetscCall(VecCopy(evec[i], vt2));
236 tmp = -eval[i];
237 PetscCall(VecAXPY(vt1, tmp, vt2));
238 PetscCall(VecNorm(vt1, NORM_INFINITY, &norm));
239 norm = PetscAbsScalar(norm);
240 if (norm > norm_max) norm_max = norm;
241 /* sniff, and bark if necessary */
242 if (norm > tols[0]) PetscCall(PetscPrintf(PETSC_COMM_SELF, " residual violation: %" PetscInt_FMT ", resi: %g\n", i, (double)norm));
243 }
244 PetscCall(PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max));
245 break;
246 default:
247 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%" PetscInt_FMT " is not supported \n", cklvl));
248 }
249 PetscCall(VecDestroy(&vt2));
250 PetscCall(VecDestroy(&vt1));
251 PetscFunctionReturn(PETSC_SUCCESS);
252 }
253
254 /*TEST
255
256 build:
257 requires: !complex
258
259 test:
260 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
261 args: -f ${DATAFILESPATH}/matrices/small
262 output_file: output/ex116_1.out
263
264 test:
265 suffix: 2
266 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
267 args: -f ${DATAFILESPATH}/matrices/small -test_syev -check_symmetry
268
269 TEST*/
270