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; 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++) { 92 PetscCall(VecResetArray(evecs[i])); 93 } 94 95 /* free space */ 96 97 PetscCall(MatDestroy(&T)); 98 99 for (i=0; i<nevs; i++) PetscCall(VecDestroy(&evecs[i])); 100 PetscCall(PetscFree(evecs)); 101 PetscCall(PetscFree(D)); 102 PetscCall(PetscFree(work)); 103 PetscCall(PetscFree(iwork)); 104 PetscCall(PetscFree(iblock)); 105 PetscCall(PetscFree(evecs_array)); 106 PetscCall(PetscFree(ifail)); 107 PetscCall(PetscFinalize()); 108 return 0; 109 #endif 110 } 111 /*------------------------------------------------ 112 Check the accuracy of the eigen solution 113 ----------------------------------------------- */ 114 /* 115 input: 116 cklvl - check level: 117 1: check residual 118 2: 1 and check B-orthogonality locally 119 A - matrix 120 il,iu - lower and upper index bound of eigenvalues 121 eval, evec - eigenvalues and eigenvectors stored in this process 122 tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] || 123 tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij 124 */ 125 #undef DEBUG_CkEigenSolutions 126 PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscScalar *eval,Vec *evec,PetscReal *tols) 127 { 128 PetscInt ierr,i,j,nev; 129 Vec vt1,vt2; /* tmp vectors */ 130 PetscReal norm,norm_max; 131 PetscScalar dot,tmp; 132 PetscReal dot_max; 133 134 PetscFunctionBegin; 135 nev = iu - il; 136 if (nev <= 0) PetscFunctionReturn(0); 137 138 PetscCall(VecDuplicate(evec[0],&vt1)); 139 PetscCall(VecDuplicate(evec[0],&vt2)); 140 141 switch (cklvl) { 142 case 2: 143 dot_max = 0.0; 144 for (i = il; i<iu; i++) { 145 PetscCall(VecCopy(evec[i], vt1)); 146 for (j=il; j<iu; j++) { 147 PetscCall(VecDot(evec[j],vt1,&dot)); 148 if (j == i) { 149 dot = PetscAbsScalar(dot - (PetscScalar)1.0); 150 } else { 151 dot = PetscAbsScalar(dot); 152 } 153 if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot); 154 #if defined(DEBUG_CkEigenSolutions) 155 if (dot > tols[1]) { 156 PetscCall(VecNorm(evec[i],NORM_INFINITY,&norm)); 157 PetscCall(PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)dot,(double)norm)); 158 } 159 #endif 160 } 161 } 162 PetscCall(PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max)); 163 164 case 1: 165 norm_max = 0.0; 166 for (i = il; i< iu; i++) { 167 PetscCall(MatMult(A, evec[i], vt1)); 168 PetscCall(VecCopy(evec[i], vt2)); 169 tmp = -eval[i]; 170 PetscCall(VecAXPY(vt1,tmp,vt2)); 171 PetscCall(VecNorm(vt1, NORM_INFINITY, &norm)); 172 norm = PetscAbsReal(norm); 173 if (norm > norm_max) norm_max = norm; 174 #if defined(DEBUG_CkEigenSolutions) 175 if (norm > tols[0]) { 176 PetscCall(PetscPrintf(PETSC_COMM_SELF," residual violation: %d, resi: %g\n",i, norm)); 177 } 178 #endif 179 } 180 PetscCall(PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max)); 181 break; 182 default: 183 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl)); 184 } 185 186 PetscCall(VecDestroy(&vt2)); 187 PetscCall(VecDestroy(&vt1)); 188 PetscFunctionReturn(0); 189 } 190