Actual source code: ex116.c
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";
10: #include <petscmat.h>
11: #include <petscblaslapack.h>
13: extern PetscErrorCode CkEigenSolutions(PetscInt, Mat, PetscInt, PetscInt, PetscReal *, Vec *, PetscReal *);
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];
32: PetscInitialize(&argc, &args, (char *)0, help);
33: MPI_Comm_size(PETSC_COMM_WORLD, &size);
36: PetscOptionsHasName(NULL, NULL, "-test_syev", &flg);
37: if (flg) TestSYEVX = PETSC_FALSE;
39: /* Determine files from which we read the two matrices */
40: PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg);
42: /* Load matrix A */
43: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[0], FILE_MODE_READ, &fd);
44: MatCreate(PETSC_COMM_WORLD, &A);
45: MatSetType(A, MATSEQAIJ);
46: MatLoad(A, fd);
47: PetscViewerDestroy(&fd);
48: MatGetSize(A, &m, &n);
50: /* Check whether A is symmetric */
51: PetscOptionsHasName(NULL, NULL, "-check_symmetry", &flg);
52: if (flg) {
53: Mat Trans;
54: MatTranspose(A, MAT_INITIAL_MATRIX, &Trans);
55: MatEqual(A, Trans, &isSymmetric);
57: MatDestroy(&Trans);
58: }
60: /* Solve eigenvalue problem: A_dense*x = lambda*B*x */
61: /*==================================================*/
62: /* Convert aij matrix to MatSeqDense for LAPACK */
63: MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense);
65: PetscBLASIntCast(8 * n, &lwork);
66: PetscBLASIntCast(n, &bn);
67: PetscMalloc1(n, &evals);
68: PetscMalloc1(lwork, &work);
69: MatDenseGetArray(A_dense, &arrayA);
71: if (!TestSYEVX) { /* test syev() */
72: 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: PetscBLASIntCast(m, &nevs);
76: il = 1;
77: PetscBLASIntCast(m, &iu);
78: } else { /* test syevx() */
79: il = 1;
80: PetscBLASIntCast(0.2 * m, &iu);
81: PetscBLASIntCast(n, &in);
82: PetscPrintf(PETSC_COMM_SELF, " LAPACKsyevx: compute %" PetscBLASInt_FMT " to %" PetscBLASInt_FMT "-th eigensolutions...\n", il, iu);
83: PetscMalloc1(m * n + 1, &evecs_array);
84: PetscMalloc1(6 * n + 1, &iwork);
85: ifail = iwork + 5 * n;
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: PetscFree(iwork);
92: }
93: MatDenseRestoreArray(A_dense, &arrayA);
96: /* View eigenvalues */
97: PetscOptionsHasName(NULL, NULL, "-eig_view", &flg);
98: if (flg) {
99: PetscPrintf(PETSC_COMM_SELF, " %" PetscBLASInt_FMT " evals: \n", nevs);
100: for (i = 0; i < nevs; i++) PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %g\n", (PetscInt)(i + il), (double)evals[i]);
101: }
103: /* Check residuals and orthogonality */
104: PetscMalloc1(nevs + 1, &evecs);
105: for (i = 0; i < nevs; i++) {
106: VecCreate(PETSC_COMM_SELF, &evecs[i]);
107: VecSetSizes(evecs[i], PETSC_DECIDE, n);
108: VecSetFromOptions(evecs[i]);
109: VecPlaceArray(evecs[i], evecs_array + i * n);
110: }
112: tols[0] = tols[1] = PETSC_SQRT_MACHINE_EPSILON;
113: CkEigenSolutions(cklvl, A, il - 1, iu - 1, evals, evecs, tols);
115: /* Free work space. */
116: for (i = 0; i < nevs; i++) VecDestroy(&evecs[i]);
117: PetscFree(evecs);
118: MatDestroy(&A_dense);
119: PetscFree(work);
120: if (TestSYEVX) PetscFree(evecs_array);
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;
133: MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense);
135: minMN = PetscMin(m, n);
136: maxMN = PetscMax(m, n);
137: lwork = 5 * minMN + maxMN;
138: PetscMalloc4(m * minMN, &arrayU, m * minMN, &arrayVT, m * minMN, &arrayErr, lwork, &work);
140: /* Create matrix Err for checking error */
141: MatCreate(PETSC_COMM_WORLD, &Err);
142: MatSetSizes(Err, PETSC_DECIDE, PETSC_DECIDE, m, minMN);
143: MatSetType(Err, MATSEQDENSE);
144: MatSeqDenseSetPreallocation(Err, (PetscScalar *)arrayErr);
146: /* Save A to arrayErr for checking accuracy later. arrayA will be destroyed by LAPACKgesvd_() */
147: MatDenseGetArray(A_dense, &arrayA);
148: PetscArraycpy(arrayErr, arrayA, m * minMN);
150: PetscBLASIntCast(m, &im);
151: 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: MatDenseRestoreArray(A_dense, &arrayA);
155: if (!lierr) {
156: PetscPrintf(PETSC_COMM_SELF, " 1st 10 of %" PetscBLASInt_FMT " singular values: \n", minMN);
157: for (i = 0; i < 10; i++) PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %g\n", i, (double)evals[i]);
158: } else {
159: PetscPrintf(PETSC_COMM_SELF, "LAPACKgesvd_ fails!");
160: }
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: MatNorm(Err, NORM_FROBENIUS, &norm);
170: PetscPrintf(PETSC_COMM_SELF, " || U*Sigma*VT - A || = %g\n", (double)norm);
172: PetscFree4(arrayU, arrayVT, arrayErr, work);
173: PetscFree(evals);
174: MatDestroy(&A_dense);
175: MatDestroy(&Err);
176: }
178: MatDestroy(&A);
179: 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: */
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;
202: nev = iu - il;
203: if (nev <= 0) return 0;
205: /*VecView(evec[0],PETSC_VIEWER_STDOUT_WORLD);*/
206: VecDuplicate(evec[0], &vt1);
207: VecDuplicate(evec[0], &vt2);
209: switch (cklvl) {
210: case 2:
211: dot_max = 0.0;
212: for (i = il; i < iu; i++) {
213: VecCopy(evec[i], vt1);
214: for (j = il; j < iu; j++) {
215: VecDot(evec[j], vt1, &dot);
216: if (j == i) {
217: dot = PetscAbsScalar(dot - 1);
218: } else {
219: dot = PetscAbsScalar(dot);
220: }
221: if (dot > dot_max) dot_max = dot;
222: if (dot > tols[1]) {
223: VecNorm(evec[i], NORM_INFINITY, &norm);
224: PetscPrintf(PETSC_COMM_SELF, "|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n", i, j, (double)dot, (double)norm);
225: }
226: }
227: }
228: PetscPrintf(PETSC_COMM_SELF, " max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max);
230: case 1:
231: norm_max = 0.0;
232: for (i = il; i < iu; i++) {
233: MatMult(A, evec[i], vt1);
234: VecCopy(evec[i], vt2);
235: tmp = -eval[i];
236: VecAXPY(vt1, tmp, vt2);
237: VecNorm(vt1, NORM_INFINITY, &norm);
238: norm = PetscAbsScalar(norm);
239: if (norm > norm_max) norm_max = norm;
240: /* sniff, and bark if necessary */
241: if (norm > tols[0]) PetscPrintf(PETSC_COMM_SELF, " residual violation: %" PetscInt_FMT ", resi: %g\n", i, (double)norm);
242: }
243: PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max);
244: break;
245: default:
246: PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%" PetscInt_FMT " is not supported \n", cklvl);
247: }
248: VecDestroy(&vt2);
249: VecDestroy(&vt1);
250: return 0;
251: }
253: /*TEST
255: build:
256: requires: !complex
258: test:
259: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
260: args: -f ${DATAFILESPATH}/matrices/small
261: output_file: output/ex116_1.out
263: test:
264: suffix: 2
265: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
266: args: -f ${DATAFILESPATH}/matrices/small -test_syev -check_symmetry
268: TEST*/