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