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 { 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, (char *)0, 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 minMN = PetscMin(m, n); 136 maxMN = PetscMax(m, n); 137 lwork = 5 * minMN + maxMN; 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 */ 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