1 static char help[] = "Test LAPACK routine DSTEBZ() and DTEIN(). \n\n"; 2 3 #include <petscmat.h> 4 #include <petscblaslapack.h> 5 6 extern PetscErrorCode CkEigenSolutions(PetscInt, Mat, PetscInt, PetscInt, PetscScalar *, Vec *, PetscReal *); 7 8 int main(int argc, char **args) { 9 #if defined(PETSC_USE_COMPLEX) || defined(PETSC_MISSING_LAPACK_STEBZ) || defined(PETSC_MISSING_LAPACK_STEIN) 10 PetscFunctionBeginUser; 11 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 12 SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP_SYS, "This example requires LAPACK routines dstebz and stien and real numbers"); 13 #else 14 PetscReal *work, tols[2]; 15 PetscInt i, j; 16 PetscBLASInt n, il = 1, iu = 5, *iblock, *isplit, *iwork, nevs, *ifail, cklvl = 2; 17 PetscMPIInt size; 18 PetscBool flg; 19 Vec *evecs; 20 PetscScalar *evecs_array, *D, *E, *evals; 21 Mat T; 22 PetscReal vl = 0.0, vu = 4.0, tol = 1000 * PETSC_MACHINE_EPSILON; 23 PetscBLASInt nsplit, info; 24 25 PetscFunctionBeginUser; 26 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 27 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 28 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 29 30 n = 100; 31 nevs = iu - il; 32 PetscCall(PetscMalloc1(3 * n + 1, &D)); 33 E = D + n; 34 evals = E + n; 35 PetscCall(PetscMalloc1(5 * n + 1, &work)); 36 PetscCall(PetscMalloc1(3 * n + 1, &iwork)); 37 PetscCall(PetscMalloc1(3 * n + 1, &iblock)); 38 isplit = iblock + n; 39 40 /* Set symmetric tridiagonal matrix */ 41 for (i = 0; i < n; i++) { 42 D[i] = 2.0; 43 E[i] = 1.0; 44 } 45 46 /* Solve eigenvalue problem: A*evec = eval*evec */ 47 PetscCall(PetscPrintf(PETSC_COMM_SELF, " LAPACKstebz_: compute %d eigenvalues...\n", nevs)); 48 LAPACKstebz_("I", "E", &n, &vl, &vu, &il, &iu, &tol, (PetscReal *)D, (PetscReal *)E, &nevs, &nsplit, (PetscReal *)evals, iblock, isplit, work, iwork, &info); 49 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_USER, "LAPACKstebz_ fails. info %d", info); 50 51 PetscCall(PetscPrintf(PETSC_COMM_SELF, " LAPACKstein_: compute %d found eigenvectors...\n", nevs)); 52 PetscCall(PetscMalloc1(n * nevs, &evecs_array)); 53 PetscCall(PetscMalloc1(nevs, &ifail)); 54 LAPACKstein_(&n, (PetscReal *)D, (PetscReal *)E, &nevs, (PetscReal *)evals, iblock, isplit, evecs_array, &n, work, iwork, ifail, &info); 55 PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_USER, "LAPACKstein_ fails. info %d", info); 56 /* View evals */ 57 PetscCall(PetscOptionsHasName(NULL, NULL, "-eig_view", &flg)); 58 if (flg) { 59 PetscCall(PetscPrintf(PETSC_COMM_SELF, " %d evals: \n", nevs)); 60 for (i = 0; i < nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT " %g\n", i, (double)evals[i])); 61 } 62 63 /* Check residuals and orthogonality */ 64 PetscCall(MatCreate(PETSC_COMM_SELF, &T)); 65 PetscCall(MatSetSizes(T, PETSC_DECIDE, PETSC_DECIDE, n, n)); 66 PetscCall(MatSetType(T, MATSBAIJ)); 67 PetscCall(MatSetFromOptions(T)); 68 PetscCall(MatSetUp(T)); 69 for (i = 0; i < n; i++) { 70 PetscCall(MatSetValues(T, 1, &i, 1, &i, &D[i], INSERT_VALUES)); 71 if (i != n - 1) { 72 j = i + 1; 73 PetscCall(MatSetValues(T, 1, &i, 1, &j, &E[i], INSERT_VALUES)); 74 } 75 } 76 PetscCall(MatAssemblyBegin(T, MAT_FINAL_ASSEMBLY)); 77 PetscCall(MatAssemblyEnd(T, MAT_FINAL_ASSEMBLY)); 78 79 PetscCall(PetscMalloc1(nevs + 1, &evecs)); 80 for (i = 0; i < nevs; i++) { 81 PetscCall(VecCreate(PETSC_COMM_SELF, &evecs[i])); 82 PetscCall(VecSetSizes(evecs[i], PETSC_DECIDE, n)); 83 PetscCall(VecSetFromOptions(evecs[i])); 84 PetscCall(VecPlaceArray(evecs[i], evecs_array + i * n)); 85 } 86 87 tols[0] = 1.e-8; 88 tols[1] = 1.e-8; 89 PetscCall(CkEigenSolutions(cklvl, T, il - 1, iu - 1, evals, evecs, tols)); 90 91 for (i = 0; i < nevs; i++) { PetscCall(VecResetArray(evecs[i])); } 92 93 /* free space */ 94 95 PetscCall(MatDestroy(&T)); 96 97 for (i = 0; i < nevs; i++) PetscCall(VecDestroy(&evecs[i])); 98 PetscCall(PetscFree(evecs)); 99 PetscCall(PetscFree(D)); 100 PetscCall(PetscFree(work)); 101 PetscCall(PetscFree(iwork)); 102 PetscCall(PetscFree(iblock)); 103 PetscCall(PetscFree(evecs_array)); 104 PetscCall(PetscFree(ifail)); 105 PetscCall(PetscFinalize()); 106 return 0; 107 #endif 108 } 109 /*------------------------------------------------ 110 Check the accuracy of the eigen solution 111 ----------------------------------------------- */ 112 /* 113 input: 114 cklvl - check level: 115 1: check residual 116 2: 1 and check B-orthogonality locally 117 A - matrix 118 il,iu - lower and upper index bound of eigenvalues 119 eval, evec - eigenvalues and eigenvectors stored in this process 120 tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] || 121 tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij 122 */ 123 #undef DEBUG_CkEigenSolutions 124 PetscErrorCode CkEigenSolutions(PetscInt cklvl, Mat A, PetscInt il, PetscInt iu, PetscScalar *eval, Vec *evec, PetscReal *tols) { 125 PetscInt ierr, i, j, nev; 126 Vec vt1, vt2; /* tmp vectors */ 127 PetscReal norm, norm_max; 128 PetscScalar dot, tmp; 129 PetscReal dot_max; 130 131 PetscFunctionBegin; 132 nev = iu - il; 133 if (nev <= 0) PetscFunctionReturn(0); 134 135 PetscCall(VecDuplicate(evec[0], &vt1)); 136 PetscCall(VecDuplicate(evec[0], &vt2)); 137 138 switch (cklvl) { 139 case 2: 140 dot_max = 0.0; 141 for (i = il; i < iu; i++) { 142 PetscCall(VecCopy(evec[i], vt1)); 143 for (j = il; j < iu; j++) { 144 PetscCall(VecDot(evec[j], vt1, &dot)); 145 if (j == i) { 146 dot = PetscAbsScalar(dot - (PetscScalar)1.0); 147 } else { 148 dot = PetscAbsScalar(dot); 149 } 150 if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot); 151 #if defined(DEBUG_CkEigenSolutions) 152 if (dot > tols[1]) { 153 PetscCall(VecNorm(evec[i], NORM_INFINITY, &norm)); 154 PetscCall(PetscPrintf(PETSC_COMM_SELF, "|delta(%d,%d)|: %g, norm: %d\n", i, j, (double)dot, (double)norm)); 155 } 156 #endif 157 } 158 } 159 PetscCall(PetscPrintf(PETSC_COMM_SELF, " max|(x_j^T*x_i) - delta_ji|: %g\n", (double)dot_max)); 160 161 case 1: 162 norm_max = 0.0; 163 for (i = il; i < iu; i++) { 164 PetscCall(MatMult(A, evec[i], vt1)); 165 PetscCall(VecCopy(evec[i], vt2)); 166 tmp = -eval[i]; 167 PetscCall(VecAXPY(vt1, tmp, vt2)); 168 PetscCall(VecNorm(vt1, NORM_INFINITY, &norm)); 169 norm = PetscAbsReal(norm); 170 if (norm > norm_max) norm_max = norm; 171 #if defined(DEBUG_CkEigenSolutions) 172 if (norm > tols[0]) { PetscCall(PetscPrintf(PETSC_COMM_SELF, " residual violation: %d, resi: %g\n", i, norm)); } 173 #endif 174 } 175 PetscCall(PetscPrintf(PETSC_COMM_SELF, " max_resi: %g\n", (double)norm_max)); 176 break; 177 default: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Error: cklvl=%d is not supported \n", cklvl)); 178 } 179 180 PetscCall(VecDestroy(&vt2)); 181 PetscCall(VecDestroy(&vt1)); 182 PetscFunctionReturn(0); 183 } 184