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 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 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 26 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 27 PetscCheck(size == 1,PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 28 29 n = 100; 30 nevs = iu - il; 31 PetscCall(PetscMalloc1(3*n+1,&D)); 32 E = D + n; 33 evals = E + n; 34 PetscCall(PetscMalloc1(5*n+1,&work)); 35 PetscCall(PetscMalloc1(3*n+1,&iwork)); 36 PetscCall(PetscMalloc1(3*n+1,&iblock)); 37 isplit = iblock + n; 38 39 /* Set symmetric tridiagonal matrix */ 40 for (i=0; i<n; i++) { 41 D[i] = 2.0; 42 E[i] = 1.0; 43 } 44 45 /* Solve eigenvalue problem: A*evec = eval*evec */ 46 PetscCall(PetscPrintf(PETSC_COMM_SELF," LAPACKstebz_: compute %d eigenvalues...\n",nevs)); 47 LAPACKstebz_("I","E",&n,&vl,&vu,&il,&iu,&tol,(PetscReal*)D,(PetscReal*)E,&nevs,&nsplit,(PetscReal*)evals,iblock,isplit,work,iwork,&info); 48 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstebz_ fails. info %d",info); 49 50 PetscCall(PetscPrintf(PETSC_COMM_SELF," LAPACKstein_: compute %d found eigenvectors...\n",nevs)); 51 PetscCall(PetscMalloc1(n*nevs,&evecs_array)); 52 PetscCall(PetscMalloc1(nevs,&ifail)); 53 LAPACKstein_(&n,(PetscReal*)D,(PetscReal*)E,&nevs,(PetscReal*)evals,iblock,isplit,evecs_array,&n,work,iwork,ifail,&info); 54 PetscCheck(!info,PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstein_ fails. info %d",info); 55 /* View evals */ 56 PetscCall(PetscOptionsHasName(NULL,NULL, "-eig_view", &flg)); 57 if (flg) { 58 PetscCall(PetscPrintf(PETSC_COMM_SELF," %d evals: \n",nevs)); 59 for (i=0; i<nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_SELF,"%" PetscInt_FMT " %g\n",i,(double)evals[i])); 60 } 61 62 /* Check residuals and orthogonality */ 63 PetscCall(MatCreate(PETSC_COMM_SELF,&T)); 64 PetscCall(MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n)); 65 PetscCall(MatSetType(T,MATSBAIJ)); 66 PetscCall(MatSetFromOptions(T)); 67 PetscCall(MatSetUp(T)); 68 for (i=0; i<n; i++) { 69 PetscCall(MatSetValues(T,1,&i,1,&i,&D[i],INSERT_VALUES)); 70 if (i != n-1) { 71 j = i+1; 72 PetscCall(MatSetValues(T,1,&i,1,&j,&E[i],INSERT_VALUES)); 73 } 74 } 75 PetscCall(MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY)); 76 PetscCall(MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY)); 77 78 PetscCall(PetscMalloc1(nevs+1,&evecs)); 79 for (i=0; i<nevs; i++) { 80 PetscCall(VecCreate(PETSC_COMM_SELF,&evecs[i])); 81 PetscCall(VecSetSizes(evecs[i],PETSC_DECIDE,n)); 82 PetscCall(VecSetFromOptions(evecs[i])); 83 PetscCall(VecPlaceArray(evecs[i],evecs_array+i*n)); 84 } 85 86 tols[0] = 1.e-8; tols[1] = 1.e-8; 87 PetscCall(CkEigenSolutions(cklvl,T,il-1,iu-1,evals,evecs,tols)); 88 89 for (i=0; i<nevs; i++) { 90 PetscCall(VecResetArray(evecs[i])); 91 } 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 { 126 PetscInt ierr,i,j,nev; 127 Vec vt1,vt2; /* tmp vectors */ 128 PetscReal norm,norm_max; 129 PetscScalar dot,tmp; 130 PetscReal dot_max; 131 132 PetscFunctionBegin; 133 nev = iu - il; 134 if (nev <= 0) PetscFunctionReturn(0); 135 136 PetscCall(VecDuplicate(evec[0],&vt1)); 137 PetscCall(VecDuplicate(evec[0],&vt2)); 138 139 switch (cklvl) { 140 case 2: 141 dot_max = 0.0; 142 for (i = il; i<iu; i++) { 143 PetscCall(VecCopy(evec[i], vt1)); 144 for (j=il; j<iu; j++) { 145 PetscCall(VecDot(evec[j],vt1,&dot)); 146 if (j == i) { 147 dot = PetscAbsScalar(dot - (PetscScalar)1.0); 148 } else { 149 dot = PetscAbsScalar(dot); 150 } 151 if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot); 152 #if defined(DEBUG_CkEigenSolutions) 153 if (dot > tols[1]) { 154 PetscCall(VecNorm(evec[i],NORM_INFINITY,&norm)); 155 PetscCall(PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)dot,(double)norm)); 156 } 157 #endif 158 } 159 } 160 PetscCall(PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max)); 161 162 case 1: 163 norm_max = 0.0; 164 for (i = il; i< iu; i++) { 165 PetscCall(MatMult(A, evec[i], vt1)); 166 PetscCall(VecCopy(evec[i], vt2)); 167 tmp = -eval[i]; 168 PetscCall(VecAXPY(vt1,tmp,vt2)); 169 PetscCall(VecNorm(vt1, NORM_INFINITY, &norm)); 170 norm = PetscAbsReal(norm); 171 if (norm > norm_max) norm_max = norm; 172 #if defined(DEBUG_CkEigenSolutions) 173 if (norm > tols[0]) { 174 PetscCall(PetscPrintf(PETSC_COMM_SELF," residual violation: %d, resi: %g\n",i, norm)); 175 } 176 #endif 177 } 178 PetscCall(PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max)); 179 break; 180 default: 181 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl)); 182 } 183 184 PetscCall(VecDestroy(&vt2)); 185 PetscCall(VecDestroy(&vt1)); 186 PetscFunctionReturn(0); 187 } 188