1 static char help[] = "Test LAPACK routine ZHEEV, ZHEEVX, ZHEGV and ZHEGVX. \n\ 2 ZHEEV computes all eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A. \n\n"; 3 4 #include <petscmat.h> 5 #include <petscblaslapack.h> 6 7 extern PetscErrorCode CkEigenSolutions(PetscInt, Mat, PetscInt, PetscInt, PetscReal *, Vec *, PetscReal *); 8 9 int main(int argc, char **args) 10 { 11 Mat A, A_dense, B; 12 Vec *evecs; 13 PetscBool flg, TestZHEEV = PETSC_TRUE, TestZHEEVX = PETSC_FALSE, TestZHEGV = PETSC_FALSE, TestZHEGVX = PETSC_FALSE; 14 PetscBool isSymmetric; 15 PetscScalar *arrayA, *arrayB, *evecs_array = NULL, *work; 16 PetscReal *evals, *rwork; 17 PetscMPIInt size; 18 PetscInt m, i, j, cklvl = 2; 19 PetscReal vl, vu, abstol = 1.e-8; 20 PetscBLASInt nn, nevs, il, iu, *iwork, *ifail, lwork, lierr, bn, one = 1; 21 PetscReal tols[2]; 22 PetscScalar v, sigma2; 23 PetscRandom rctx; 24 PetscReal h2, sigma1 = 100.0; 25 PetscInt dim, Ii, J, n = 6, use_random; 26 27 PetscFunctionBeginUser; 28 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 29 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 30 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 31 32 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_zheevx", &flg)); 33 if (flg) { 34 TestZHEEV = PETSC_FALSE; 35 TestZHEEVX = PETSC_TRUE; 36 } 37 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_zhegv", &flg)); 38 if (flg) { 39 TestZHEEV = PETSC_FALSE; 40 TestZHEGV = PETSC_TRUE; 41 } 42 PetscCall(PetscOptionsHasName(NULL, NULL, "-test_zhegvx", &flg)); 43 if (flg) { 44 TestZHEEV = PETSC_FALSE; 45 TestZHEGVX = PETSC_TRUE; 46 } 47 48 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma1", &sigma1, NULL)); 49 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 50 dim = n * n; 51 52 PetscCall(MatCreate(PETSC_COMM_SELF, &A)); 53 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim)); 54 PetscCall(MatSetType(A, MATSEQDENSE)); 55 PetscCall(MatSetFromOptions(A)); 56 PetscCall(MatSetUp(A)); 57 58 PetscCall(PetscOptionsHasName(NULL, NULL, "-norandom", &flg)); 59 if (flg) use_random = 0; 60 else use_random = 1; 61 if (use_random) { 62 PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rctx)); 63 PetscCall(PetscRandomSetFromOptions(rctx)); 64 PetscCall(PetscRandomSetInterval(rctx, 0.0, PETSC_i)); 65 } else { 66 sigma2 = 10.0 * PETSC_i; 67 } 68 h2 = 1.0 / ((n + 1) * (n + 1)); 69 for (Ii = 0; Ii < dim; Ii++) { 70 v = -1.0; 71 i = Ii / n; 72 j = Ii - i * n; 73 if (i > 0) { 74 J = Ii - n; 75 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 76 } 77 if (i < n - 1) { 78 J = Ii + n; 79 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 80 } 81 if (j > 0) { 82 J = Ii - 1; 83 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 84 } 85 if (j < n - 1) { 86 J = Ii + 1; 87 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 88 } 89 if (use_random) PetscCall(PetscRandomGetValue(rctx, &sigma2)); 90 v = 4.0 - sigma1 * h2; 91 PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, ADD_VALUES)); 92 } 93 /* make A complex Hermitian */ 94 v = sigma2 * h2; 95 Ii = 0; 96 J = 1; 97 PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, ADD_VALUES)); 98 v = -sigma2 * h2; 99 PetscCall(MatSetValues(A, 1, &J, 1, &Ii, &v, ADD_VALUES)); 100 if (use_random) PetscCall(PetscRandomDestroy(&rctx)); 101 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 102 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 103 m = n = dim; 104 105 /* Check whether A is symmetric */ 106 PetscCall(PetscOptionsHasName(NULL, NULL, "-check_symmetry", &flg)); 107 if (flg) { 108 Mat Trans; 109 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Trans)); 110 PetscCall(MatEqual(A, Trans, &isSymmetric)); 111 PetscCheck(isSymmetric, PETSC_COMM_SELF, PETSC_ERR_USER, "A must be symmetric"); 112 PetscCall(MatDestroy(&Trans)); 113 } 114 115 /* Convert aij matrix to MatSeqDense for LAPACK */ 116 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQDENSE, &flg)); 117 if (flg) { 118 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &A_dense)); 119 } else { 120 PetscCall(MatConvert(A, MATSEQDENSE, MAT_INITIAL_MATRIX, &A_dense)); 121 } 122 123 PetscCall(MatCreate(PETSC_COMM_SELF, &B)); 124 PetscCall(MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, dim, dim)); 125 PetscCall(MatSetType(B, MATSEQDENSE)); 126 PetscCall(MatSetFromOptions(B)); 127 PetscCall(MatSetUp(B)); 128 v = 1.0; 129 for (Ii = 0; Ii < dim; Ii++) PetscCall(MatSetValues(B, 1, &Ii, 1, &Ii, &v, ADD_VALUES)); 130 131 /* Solve standard eigenvalue problem: A*x = lambda*x */ 132 /*===================================================*/ 133 PetscCall(PetscBLASIntCast(2 * n, &lwork)); 134 PetscCall(PetscBLASIntCast(n, &bn)); 135 PetscCall(PetscMalloc1(n, &evals)); 136 PetscCall(PetscMalloc1(lwork, &work)); 137 PetscCall(MatDenseGetArray(A_dense, &arrayA)); 138 139 if (TestZHEEV) { /* test zheev() */ 140 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n", m)); 141 PetscCall(PetscMalloc1(3 * n - 2, &rwork)); 142 LAPACKsyev_("V", "U", &bn, arrayA, &bn, evals, work, &lwork, rwork, &lierr); 143 PetscCall(PetscFree(rwork)); 144 145 evecs_array = arrayA; 146 nevs = m; 147 il = 1; 148 iu = m; 149 } 150 if (TestZHEEVX) { 151 il = 1; 152 PetscCall(PetscBLASIntCast((0.2 * m), &iu)); 153 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsyevx: compute %d to %d-th eigensolutions...\n", il, iu)); 154 PetscCall(PetscMalloc1(m * n + 1, &evecs_array)); 155 PetscCall(PetscMalloc1(7 * n + 1, &rwork)); 156 PetscCall(PetscMalloc1(5 * n + 1, &iwork)); 157 PetscCall(PetscMalloc1(n + 1, &ifail)); 158 159 /* in the case "I", vl and vu are not referenced */ 160 vl = 0.0; 161 vu = 8.0; 162 PetscCall(PetscBLASIntCast(n, &nn)); 163 LAPACKsyevx_("V", "I", "U", &bn, arrayA, &bn, &vl, &vu, &il, &iu, &abstol, &nevs, evals, evecs_array, &nn, work, &lwork, rwork, iwork, ifail, &lierr); 164 PetscCall(PetscFree(iwork)); 165 PetscCall(PetscFree(ifail)); 166 PetscCall(PetscFree(rwork)); 167 } 168 if (TestZHEGV) { 169 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsygv: compute all %" PetscInt_FMT " eigensolutions...\n", m)); 170 PetscCall(PetscMalloc1(3 * n + 1, &rwork)); 171 PetscCall(MatDenseGetArray(B, &arrayB)); 172 LAPACKsygv_(&one, "V", "U", &bn, arrayA, &bn, arrayB, &bn, evals, work, &lwork, rwork, &lierr); 173 evecs_array = arrayA; 174 nevs = m; 175 il = 1; 176 iu = m; 177 PetscCall(MatDenseRestoreArray(B, &arrayB)); 178 PetscCall(PetscFree(rwork)); 179 } 180 if (TestZHEGVX) { 181 il = 1; 182 PetscCall(PetscBLASIntCast((0.2 * m), &iu)); 183 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " LAPACKsygv: compute %d to %d-th eigensolutions...\n", il, iu)); 184 PetscCall(PetscMalloc1(m * n + 1, &evecs_array)); 185 PetscCall(PetscMalloc1(6 * n + 1, &iwork)); 186 ifail = iwork + 5 * n; 187 PetscCall(PetscMalloc1(7 * n + 1, &rwork)); 188 PetscCall(MatDenseGetArray(B, &arrayB)); 189 vl = 0.0; 190 vu = 8.0; 191 PetscCall(PetscBLASIntCast(n, &nn)); 192 LAPACKsygvx_(&one, "V", "I", "U", &bn, arrayA, &bn, arrayB, &bn, &vl, &vu, &il, &iu, &abstol, &nevs, evals, evecs_array, &nn, work, &lwork, rwork, iwork, ifail, &lierr); 193 PetscCall(MatDenseRestoreArray(B, &arrayB)); 194 PetscCall(PetscFree(iwork)); 195 PetscCall(PetscFree(rwork)); 196 } 197 PetscCall(MatDenseRestoreArray(A_dense, &arrayA)); 198 PetscCheck(nevs > 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs); 199 200 /* View evals */ 201 PetscCall(PetscOptionsHasName(NULL, NULL, "-eig_view", &flg)); 202 if (flg) { 203 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " %d evals: \n", nevs)); 204 for (i = 0; i < nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT " %g\n", i + il, (double)evals[i])); 205 } 206 207 /* Check residuals and orthogonality */ 208 PetscCall(PetscMalloc1(nevs + 1, &evecs)); 209 for (i = 0; i < nevs; i++) { 210 PetscCall(VecCreate(PETSC_COMM_SELF, &evecs[i])); 211 PetscCall(VecSetSizes(evecs[i], PETSC_DECIDE, n)); 212 PetscCall(VecSetFromOptions(evecs[i])); 213 PetscCall(VecPlaceArray(evecs[i], evecs_array + i * n)); 214 } 215 216 tols[0] = PETSC_SQRT_MACHINE_EPSILON; 217 tols[1] = PETSC_SQRT_MACHINE_EPSILON; 218 PetscCall(CkEigenSolutions(cklvl, A, il - 1, iu - 1, evals, evecs, tols)); 219 for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i])); 220 PetscCall(PetscFree(evecs)); 221 222 /* Free work space. */ 223 if (TestZHEEVX || TestZHEGVX) PetscCall(PetscFree(evecs_array)); 224 PetscCall(PetscFree(evals)); 225 PetscCall(PetscFree(work)); 226 PetscCall(MatDestroy(&A_dense)); 227 PetscCall(MatDestroy(&A)); 228 PetscCall(MatDestroy(&B)); 229 PetscCall(PetscFinalize()); 230 return 0; 231 } 232 /*------------------------------------------------ 233 Check the accuracy of the eigen solution 234 ----------------------------------------------- */ 235 /* 236 input: 237 cklvl - check level: 238 1: check residual 239 2: 1 and check B-orthogonality locally 240 A - matrix 241 il,iu - lower and upper index bound of eigenvalues 242 eval, evec - eigenvalues and eigenvectors stored in this process 243 tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] || 244 tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij 245 */ 246 PetscErrorCode CkEigenSolutions(PetscInt cklvl, Mat A, PetscInt il, PetscInt iu, PetscReal *eval, Vec *evec, PetscReal *tols) 247 { 248 PetscInt i, j, nev; 249 Vec vt1, vt2; /* tmp vectors */ 250 PetscReal norm, tmp, norm_max, dot_max, rdot; 251 PetscScalar dot; 252 253 PetscFunctionBegin; 254 nev = iu - il; 255 if (nev <= 0) PetscFunctionReturn(0); 256 257 PetscCall(VecDuplicate(evec[0], &vt1)); 258 PetscCall(VecDuplicate(evec[0], &vt2)); 259 260 switch (cklvl) { 261 case 2: 262 dot_max = 0.0; 263 for (i = il; i < iu; i++) { 264 PetscCall(VecCopy(evec[i], vt1)); 265 for (j = il; j < iu; j++) { 266 PetscCall(VecDot(evec[j], vt1, &dot)); 267 if (j == i) { 268 rdot = PetscAbsScalar(dot - (PetscScalar)1.0); 269 } else { 270 rdot = PetscAbsScalar(dot); 271 } 272 if (rdot > dot_max) dot_max = rdot; 273 if (rdot > tols[1]) { 274 PetscCall(VecNorm(evec[i], NORM_INFINITY, &norm)); 275 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n", i, j, (double)rdot, (double)norm)); 276 } 277 } 278 } 279 PetscCall(PetscPrintf(PETSC_COMM_SELF, " max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max)); 280 281 case 1: 282 norm_max = 0.0; 283 for (i = il; i < iu; i++) { 284 PetscCall(MatMult(A, evec[i], vt1)); 285 PetscCall(VecCopy(evec[i], vt2)); 286 tmp = -eval[i]; 287 PetscCall(VecAXPY(vt1, tmp, vt2)); 288 PetscCall(VecNorm(vt1, NORM_INFINITY, &norm)); 289 norm = PetscAbs(norm); 290 if (norm > norm_max) norm_max = norm; 291 /* sniff, and bark if necessary */ 292 if (norm > tols[0]) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " residual violation: %" PetscInt_FMT ", resi: %g\n", i, (double)norm)); 293 } 294 PetscCall(PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max)); 295 break; 296 default: 297 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%" PetscInt_FMT " is not supported \n", cklvl)); 298 } 299 PetscCall(VecDestroy(&vt2)); 300 PetscCall(VecDestroy(&vt1)); 301 PetscFunctionReturn(0); 302 } 303 304 /*TEST 305 306 build: 307 requires: complex 308 309 test: 310 311 test: 312 suffix: 2 313 args: -test_zheevx 314 315 test: 316 suffix: 3 317 args: -test_zhegv 318 319 test: 320 suffix: 4 321 args: -test_zhegvx 322 323 TEST*/ 324