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