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 PetscErrorCode ierr; 11 #if defined(PETSC_USE_COMPLEX) || defined(PETSC_MISSING_LAPACK_STEBZ) || defined(PETSC_MISSING_LAPACK_STEIN) 12 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 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 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 27 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 28 PetscCheckFalse(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 CHKERRQ(PetscMalloc1(3*n+1,&D)); 33 E = D + n; 34 evals = E + n; 35 CHKERRQ(PetscMalloc1(5*n+1,&work)); 36 CHKERRQ(PetscMalloc1(3*n+1,&iwork)); 37 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF," LAPACKstein_: compute %d found eigenvectors...\n",nevs)); 52 CHKERRQ(PetscMalloc1(n*nevs,&evecs_array)); 53 CHKERRQ(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 CHKERRQ(PetscOptionsHasName(NULL,NULL, "-eig_view", &flg)); 58 if (flg) { 59 CHKERRQ(PetscPrintf(PETSC_COMM_SELF," %d evals: \n",nevs)); 60 for (i=0; i<nevs; i++) CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"%" PetscInt_FMT " %g\n",i,(double)evals[i])); 61 } 62 63 /* Check residuals and orthogonality */ 64 CHKERRQ(MatCreate(PETSC_COMM_SELF,&T)); 65 CHKERRQ(MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n)); 66 CHKERRQ(MatSetType(T,MATSBAIJ)); 67 CHKERRQ(MatSetFromOptions(T)); 68 CHKERRQ(MatSetUp(T)); 69 for (i=0; i<n; i++) { 70 CHKERRQ(MatSetValues(T,1,&i,1,&i,&D[i],INSERT_VALUES)); 71 if (i != n-1) { 72 j = i+1; 73 CHKERRQ(MatSetValues(T,1,&i,1,&j,&E[i],INSERT_VALUES)); 74 } 75 } 76 CHKERRQ(MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY)); 77 CHKERRQ(MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY)); 78 79 CHKERRQ(PetscMalloc1(nevs+1,&evecs)); 80 for (i=0; i<nevs; i++) { 81 CHKERRQ(VecCreate(PETSC_COMM_SELF,&evecs[i])); 82 CHKERRQ(VecSetSizes(evecs[i],PETSC_DECIDE,n)); 83 CHKERRQ(VecSetFromOptions(evecs[i])); 84 CHKERRQ(VecPlaceArray(evecs[i],evecs_array+i*n)); 85 } 86 87 tols[0] = 1.e-8; tols[1] = 1.e-8; 88 CHKERRQ(CkEigenSolutions(cklvl,T,il-1,iu-1,evals,evecs,tols)); 89 90 for (i=0; i<nevs; i++) { 91 CHKERRQ(VecResetArray(evecs[i])); 92 } 93 94 /* free space */ 95 96 CHKERRQ(MatDestroy(&T)); 97 98 for (i=0; i<nevs; i++) CHKERRQ(VecDestroy(&evecs[i])); 99 CHKERRQ(PetscFree(evecs)); 100 CHKERRQ(PetscFree(D)); 101 CHKERRQ(PetscFree(work)); 102 CHKERRQ(PetscFree(iwork)); 103 CHKERRQ(PetscFree(iblock)); 104 CHKERRQ(PetscFree(evecs_array)); 105 CHKERRQ(PetscFree(ifail)); 106 ierr = PetscFinalize(); 107 return ierr; 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(0); 136 137 CHKERRQ(VecDuplicate(evec[0],&vt1)); 138 CHKERRQ(VecDuplicate(evec[0],&vt2)); 139 140 switch (cklvl) { 141 case 2: 142 dot_max = 0.0; 143 for (i = il; i<iu; i++) { 144 CHKERRQ(VecCopy(evec[i], vt1)); 145 for (j=il; j<iu; j++) { 146 CHKERRQ(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 CHKERRQ(VecNorm(evec[i],NORM_INFINITY,&norm)); 156 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)dot,(double)norm)); 157 } 158 #endif 159 } 160 } 161 CHKERRQ(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 CHKERRQ(MatMult(A, evec[i], vt1)); 167 CHKERRQ(VecCopy(evec[i], vt2)); 168 tmp = -eval[i]; 169 CHKERRQ(VecAXPY(vt1,tmp,vt2)); 170 CHKERRQ(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]) { 175 CHKERRQ(PetscPrintf(PETSC_COMM_SELF," residual violation: %d, resi: %g\n",i, norm)); 176 } 177 #endif 178 } 179 CHKERRQ(PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max)); 180 break; 181 default: 182 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl)); 183 } 184 185 CHKERRQ(VecDestroy(&vt2)); 186 CHKERRQ(VecDestroy(&vt1)); 187 PetscFunctionReturn(0); 188 } 189