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