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 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 28 if (size != 1) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 29 30 n = 100; 31 nevs = iu - il; 32 ierr = PetscMalloc1(3*n+1,&D);CHKERRQ(ierr); 33 E = D + n; 34 evals = E + n; 35 ierr = PetscMalloc1(5*n+1,&work);CHKERRQ(ierr); 36 ierr = PetscMalloc1(3*n+1,&iwork);CHKERRQ(ierr); 37 ierr = PetscMalloc1(3*n+1,&iblock);CHKERRQ(ierr); 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 ierr = PetscPrintf(PETSC_COMM_SELF," LAPACKstebz_: compute %d eigenvalues...\n",nevs);CHKERRQ(ierr); 48 LAPACKstebz_("I","E",&n,&vl,&vu,&il,&iu,&tol,(PetscReal*)D,(PetscReal*)E,&nevs,&nsplit,(PetscReal*)evals,iblock,isplit,work,iwork,&info); 49 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstebz_ fails. info %d",info); 50 51 ierr = PetscPrintf(PETSC_COMM_SELF," LAPACKstein_: compute %d found eigenvectors...\n",nevs);CHKERRQ(ierr); 52 ierr = PetscMalloc1(n*nevs,&evecs_array);CHKERRQ(ierr); 53 ierr = PetscMalloc1(nevs,&ifail);CHKERRQ(ierr); 54 LAPACKstein_(&n,(PetscReal*)D,(PetscReal*)E,&nevs,(PetscReal*)evals,iblock,isplit,evecs_array,&n,work,iwork,ifail,&info); 55 if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstein_ fails. info %d",info); 56 /* View evals */ 57 ierr = PetscOptionsHasName(NULL,NULL, "-eig_view", &flg);CHKERRQ(ierr); 58 if (flg) { 59 ierr = PetscPrintf(PETSC_COMM_SELF," %d evals: \n",nevs);CHKERRQ(ierr); 60 for (i=0; i<nevs; i++) {ierr = PetscPrintf(PETSC_COMM_SELF,"%D %g\n",i,(double)evals[i]);CHKERRQ(ierr);} 61 } 62 63 /* Check residuals and orthogonality */ 64 ierr = MatCreate(PETSC_COMM_SELF,&T);CHKERRQ(ierr); 65 ierr = MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr); 66 ierr = MatSetType(T,MATSBAIJ);CHKERRQ(ierr); 67 ierr = MatSetFromOptions(T);CHKERRQ(ierr); 68 ierr = MatSetUp(T);CHKERRQ(ierr); 69 for (i=0; i<n; i++) { 70 ierr = MatSetValues(T,1,&i,1,&i,&D[i],INSERT_VALUES);CHKERRQ(ierr); 71 if (i != n-1) { 72 j = i+1; 73 ierr = MatSetValues(T,1,&i,1,&j,&E[i],INSERT_VALUES);CHKERRQ(ierr); 74 } 75 } 76 ierr = MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 77 ierr = MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 78 79 ierr = PetscMalloc1(nevs+1,&evecs);CHKERRQ(ierr); 80 for (i=0; i<nevs; i++) { 81 ierr = VecCreate(PETSC_COMM_SELF,&evecs[i]);CHKERRQ(ierr); 82 ierr = VecSetSizes(evecs[i],PETSC_DECIDE,n);CHKERRQ(ierr); 83 ierr = VecSetFromOptions(evecs[i]);CHKERRQ(ierr); 84 ierr = VecPlaceArray(evecs[i],evecs_array+i*n);CHKERRQ(ierr); 85 } 86 87 tols[0] = 1.e-8; tols[1] = 1.e-8; 88 ierr = CkEigenSolutions(cklvl,T,il-1,iu-1,evals,evecs,tols);CHKERRQ(ierr); 89 90 for (i=0; i<nevs; i++) { 91 ierr = VecResetArray(evecs[i]);CHKERRQ(ierr); 92 } 93 94 /* free space */ 95 96 ierr = MatDestroy(&T);CHKERRQ(ierr); 97 98 for (i=0; i<nevs; i++) { ierr = VecDestroy(&evecs[i]);CHKERRQ(ierr);} 99 ierr = PetscFree(evecs);CHKERRQ(ierr); 100 ierr = PetscFree(D);CHKERRQ(ierr); 101 ierr = PetscFree(work);CHKERRQ(ierr); 102 ierr = PetscFree(iwork);CHKERRQ(ierr); 103 ierr = PetscFree(iblock);CHKERRQ(ierr); 104 ierr = PetscFree(evecs_array);CHKERRQ(ierr); 105 ierr = PetscFree(ifail);CHKERRQ(ierr); 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 ierr = VecDuplicate(evec[0],&vt1);CHKERRQ(ierr); 138 ierr = VecDuplicate(evec[0],&vt2);CHKERRQ(ierr); 139 140 switch (cklvl) { 141 case 2: 142 dot_max = 0.0; 143 for (i = il; i<iu; i++) { 144 ierr = VecCopy(evec[i], vt1);CHKERRQ(ierr); 145 for (j=il; j<iu; j++) { 146 ierr = VecDot(evec[j],vt1,&dot);CHKERRQ(ierr); 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 ierr = VecNorm(evec[i],NORM_INFINITY,&norm);CHKERRQ(ierr); 156 ierr = PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)dot,(double)norm);CHKERRQ(ierr); 157 } 158 #endif 159 } 160 } 161 ierr = PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);CHKERRQ(ierr); 162 163 case 1: 164 norm_max = 0.0; 165 for (i = il; i< iu; i++) { 166 ierr = MatMult(A, evec[i], vt1);CHKERRQ(ierr); 167 ierr = VecCopy(evec[i], vt2);CHKERRQ(ierr); 168 tmp = -eval[i]; 169 ierr = VecAXPY(vt1,tmp,vt2);CHKERRQ(ierr); 170 ierr = VecNorm(vt1, NORM_INFINITY, &norm);CHKERRQ(ierr); 171 norm = PetscAbsReal(norm); 172 if (norm > norm_max) norm_max = norm; 173 #if defined(DEBUG_CkEigenSolutions) 174 if (norm > tols[0]) { 175 ierr = PetscPrintf(PETSC_COMM_SELF," residual violation: %d, resi: %g\n",i, norm);CHKERRQ(ierr); 176 } 177 #endif 178 } 179 ierr = PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max);CHKERRQ(ierr); 180 break; 181 default: 182 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);CHKERRQ(ierr); 183 } 184 185 ierr = VecDestroy(&vt2);CHKERRQ(ierr); 186 ierr = VecDestroy(&vt1);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189