1 static char help[] = "Test LAPACK routine ZHEEV, ZHEEVX, ZHEGV and ZHEGVX. \n\ 2 ZHEEV computes all eigenvalues and, optionally, eigenvectors of a complex Hermitian matrix A. \n\n"; 3 4 #include <petscmat.h> 5 #include <petscblaslapack.h> 6 7 extern PetscErrorCode CkEigenSolutions(PetscInt,Mat,PetscInt,PetscInt,PetscReal*,Vec*,PetscReal*); 8 9 int main(int argc,char **args) 10 { 11 Mat A,A_dense,B; 12 Vec *evecs; 13 PetscBool flg,TestZHEEV=PETSC_TRUE,TestZHEEVX=PETSC_FALSE,TestZHEGV=PETSC_FALSE,TestZHEGVX=PETSC_FALSE; 14 PetscBool isSymmetric; 15 PetscScalar *arrayA,*arrayB,*evecs_array=NULL,*work; 16 PetscReal *evals,*rwork; 17 PetscMPIInt size; 18 PetscInt m,i,j,cklvl=2; 19 PetscReal vl,vu,abstol=1.e-8; 20 PetscBLASInt nn,nevs,il,iu,*iwork,*ifail,lwork,lierr,bn,one=1; 21 PetscReal tols[2]; 22 PetscScalar v,sigma2; 23 PetscRandom rctx; 24 PetscReal h2,sigma1 = 100.0; 25 PetscInt dim,Ii,J,n = 6,use_random; 26 27 PetscFunctionBeginUser; 28 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 29 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 30 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 31 32 PetscCall(PetscOptionsHasName(NULL,NULL, "-test_zheevx", &flg)); 33 if (flg) { 34 TestZHEEV = PETSC_FALSE; 35 TestZHEEVX = PETSC_TRUE; 36 } 37 PetscCall(PetscOptionsHasName(NULL,NULL, "-test_zhegv", &flg)); 38 if (flg) { 39 TestZHEEV = PETSC_FALSE; 40 TestZHEGV = PETSC_TRUE; 41 } 42 PetscCall(PetscOptionsHasName(NULL,NULL, "-test_zhegvx", &flg)); 43 if (flg) { 44 TestZHEEV = PETSC_FALSE; 45 TestZHEGVX = PETSC_TRUE; 46 } 47 48 PetscCall(PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL)); 49 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 50 dim = n*n; 51 52 PetscCall(MatCreate(PETSC_COMM_SELF,&A)); 53 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim)); 54 PetscCall(MatSetType(A,MATSEQDENSE)); 55 PetscCall(MatSetFromOptions(A)); 56 PetscCall(MatSetUp(A)); 57 58 PetscCall(PetscOptionsHasName(NULL,NULL,"-norandom",&flg)); 59 if (flg) use_random = 0; 60 else use_random = 1; 61 if (use_random) { 62 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rctx)); 63 PetscCall(PetscRandomSetFromOptions(rctx)); 64 PetscCall(PetscRandomSetInterval(rctx,0.0,PETSC_i)); 65 } else { 66 sigma2 = 10.0*PETSC_i; 67 } 68 h2 = 1.0/((n+1)*(n+1)); 69 for (Ii=0; Ii<dim; Ii++) { 70 v = -1.0; i = Ii/n; j = Ii - i*n; 71 if (i>0) { 72 J = Ii-n; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 73 } 74 if (i<n-1) { 75 J = Ii+n; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 76 } 77 if (j>0) { 78 J = Ii-1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 79 } 80 if (j<n-1) { 81 J = Ii+1; PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 82 } 83 if (use_random) PetscCall(PetscRandomGetValue(rctx,&sigma2)); 84 v = 4.0 - sigma1*h2; 85 PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES)); 86 } 87 /* make A complex Hermitian */ 88 v = sigma2*h2; 89 Ii = 0; J = 1; 90 PetscCall(MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES)); 91 v = -sigma2*h2; 92 PetscCall(MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES)); 93 if (use_random) PetscCall(PetscRandomDestroy(&rctx)); 94 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 95 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 96 m = n = dim; 97 98 /* Check whether A is symmetric */ 99 PetscCall(PetscOptionsHasName(NULL,NULL, "-check_symmetry", &flg)); 100 if (flg) { 101 Mat Trans; 102 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX, &Trans)); 103 PetscCall(MatEqual(A, Trans, &isSymmetric)); 104 PetscCheck(isSymmetric,PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric"); 105 PetscCall(MatDestroy(&Trans)); 106 } 107 108 /* Convert aij matrix to MatSeqDense for LAPACK */ 109 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg)); 110 if (flg) { 111 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A_dense)); 112 } else { 113 PetscCall(MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense)); 114 } 115 116 PetscCall(MatCreate(PETSC_COMM_SELF,&B)); 117 PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,dim,dim)); 118 PetscCall(MatSetType(B,MATSEQDENSE)); 119 PetscCall(MatSetFromOptions(B)); 120 PetscCall(MatSetUp(B)); 121 v = 1.0; 122 for (Ii=0; Ii<dim; Ii++) { 123 PetscCall(MatSetValues(B,1,&Ii,1,&Ii,&v,ADD_VALUES)); 124 } 125 126 /* Solve standard eigenvalue problem: A*x = lambda*x */ 127 /*===================================================*/ 128 PetscCall(PetscBLASIntCast(2*n,&lwork)); 129 PetscCall(PetscBLASIntCast(n,&bn)); 130 PetscCall(PetscMalloc1(n,&evals)); 131 PetscCall(PetscMalloc1(lwork,&work)); 132 PetscCall(MatDenseGetArray(A_dense,&arrayA)); 133 134 if (TestZHEEV) { /* test zheev() */ 135 PetscCall(PetscPrintf(PETSC_COMM_WORLD," LAPACKsyev: compute all %" PetscInt_FMT " eigensolutions...\n",m)); 136 PetscCall(PetscMalloc1(3*n-2,&rwork)); 137 LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,rwork,&lierr); 138 PetscCall(PetscFree(rwork)); 139 140 evecs_array = arrayA; 141 nevs = m; 142 il =1; iu=m; 143 } 144 if (TestZHEEVX) { 145 il = 1; 146 PetscCall(PetscBLASIntCast((0.2*m),&iu)); 147 PetscCall(PetscPrintf(PETSC_COMM_WORLD," LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu)); 148 PetscCall(PetscMalloc1(m*n+1,&evecs_array)); 149 PetscCall(PetscMalloc1(7*n+1,&rwork)); 150 PetscCall(PetscMalloc1(5*n+1,&iwork)); 151 PetscCall(PetscMalloc1(n+1,&ifail)); 152 153 /* in the case "I", vl and vu are not referenced */ 154 vl = 0.0; vu = 8.0; 155 PetscCall(PetscBLASIntCast(n,&nn)); 156 LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&nn,work,&lwork,rwork,iwork,ifail,&lierr); 157 PetscCall(PetscFree(iwork)); 158 PetscCall(PetscFree(ifail)); 159 PetscCall(PetscFree(rwork)); 160 } 161 if (TestZHEGV) { 162 PetscCall(PetscPrintf(PETSC_COMM_WORLD," LAPACKsygv: compute all %" PetscInt_FMT " eigensolutions...\n",m)); 163 PetscCall(PetscMalloc1(3*n+1,&rwork)); 164 PetscCall(MatDenseGetArray(B,&arrayB)); 165 LAPACKsygv_(&one,"V","U",&bn,arrayA,&bn,arrayB,&bn,evals,work,&lwork,rwork,&lierr); 166 evecs_array = arrayA; 167 nevs = m; 168 il = 1; iu=m; 169 PetscCall(MatDenseRestoreArray(B,&arrayB)); 170 PetscCall(PetscFree(rwork)); 171 } 172 if (TestZHEGVX) { 173 il = 1; 174 PetscCall(PetscBLASIntCast((0.2*m),&iu)); 175 PetscCall(PetscPrintf(PETSC_COMM_WORLD," LAPACKsygv: compute %d to %d-th eigensolutions...\n",il,iu)); 176 PetscCall(PetscMalloc1(m*n+1,&evecs_array)); 177 PetscCall(PetscMalloc1(6*n+1,&iwork)); 178 ifail = iwork + 5*n; 179 PetscCall(PetscMalloc1(7*n+1,&rwork)); 180 PetscCall(MatDenseGetArray(B,&arrayB)); 181 vl = 0.0; vu = 8.0; 182 PetscCall(PetscBLASIntCast(n,&nn)); 183 LAPACKsygvx_(&one,"V","I","U",&bn,arrayA,&bn,arrayB,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&nn,work,&lwork,rwork,iwork,ifail,&lierr); 184 PetscCall(MatDenseRestoreArray(B,&arrayB)); 185 PetscCall(PetscFree(iwork)); 186 PetscCall(PetscFree(rwork)); 187 } 188 PetscCall(MatDenseRestoreArray(A_dense,&arrayA)); 189 PetscCheck(nevs > 0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs); 190 191 /* View evals */ 192 PetscCall(PetscOptionsHasName(NULL,NULL, "-eig_view", &flg)); 193 if (flg) { 194 PetscCall(PetscPrintf(PETSC_COMM_WORLD," %d evals: \n",nevs)); 195 for (i=0; i<nevs; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT " %g\n",i+il,(double)evals[i])); 196 } 197 198 /* Check residuals and orthogonality */ 199 PetscCall(PetscMalloc1(nevs+1,&evecs)); 200 for (i=0; i<nevs; i++) { 201 PetscCall(VecCreate(PETSC_COMM_SELF,&evecs[i])); 202 PetscCall(VecSetSizes(evecs[i],PETSC_DECIDE,n)); 203 PetscCall(VecSetFromOptions(evecs[i])); 204 PetscCall(VecPlaceArray(evecs[i],evecs_array+i*n)); 205 } 206 207 tols[0] = PETSC_SQRT_MACHINE_EPSILON; tols[1] = PETSC_SQRT_MACHINE_EPSILON; 208 PetscCall(CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols)); 209 for (i=0; i<nevs; i++) PetscCall(VecDestroy(&evecs[i])); 210 PetscCall(PetscFree(evecs)); 211 212 /* Free work space. */ 213 if (TestZHEEVX || TestZHEGVX) { 214 PetscCall(PetscFree(evecs_array)); 215 } 216 PetscCall(PetscFree(evals)); 217 PetscCall(PetscFree(work)); 218 PetscCall(MatDestroy(&A_dense)); 219 PetscCall(MatDestroy(&A)); 220 PetscCall(MatDestroy(&B)); 221 PetscCall(PetscFinalize()); 222 return 0; 223 } 224 /*------------------------------------------------ 225 Check the accuracy of the eigen solution 226 ----------------------------------------------- */ 227 /* 228 input: 229 cklvl - check level: 230 1: check residual 231 2: 1 and check B-orthogonality locally 232 A - matrix 233 il,iu - lower and upper index bound of eigenvalues 234 eval, evec - eigenvalues and eigenvectors stored in this process 235 tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] || 236 tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij 237 */ 238 PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols) 239 { 240 PetscInt i,j,nev; 241 Vec vt1,vt2; /* tmp vectors */ 242 PetscReal norm,tmp,norm_max,dot_max,rdot; 243 PetscScalar dot; 244 245 PetscFunctionBegin; 246 nev = iu - il; 247 if (nev <= 0) PetscFunctionReturn(0); 248 249 PetscCall(VecDuplicate(evec[0],&vt1)); 250 PetscCall(VecDuplicate(evec[0],&vt2)); 251 252 switch (cklvl) { 253 case 2: 254 dot_max = 0.0; 255 for (i = il; i<iu; i++) { 256 PetscCall(VecCopy(evec[i], vt1)); 257 for (j=il; j<iu; j++) { 258 PetscCall(VecDot(evec[j],vt1,&dot)); 259 if (j == i) { 260 rdot = PetscAbsScalar(dot - (PetscScalar)1.0); 261 } else { 262 rdot = PetscAbsScalar(dot); 263 } 264 if (rdot > dot_max) dot_max = rdot; 265 if (rdot > tols[1]) { 266 PetscCall(VecNorm(evec[i],NORM_INFINITY,&norm)); 267 PetscCall(PetscPrintf(PETSC_COMM_SELF,"|delta(%" PetscInt_FMT ",%" PetscInt_FMT ")|: %g, norm: %g\n",i,j,(double)rdot,(double)norm)); 268 } 269 } 270 } 271 PetscCall(PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max)); 272 273 case 1: 274 norm_max = 0.0; 275 for (i = il; i< iu; i++) { 276 PetscCall(MatMult(A, evec[i], vt1)); 277 PetscCall(VecCopy(evec[i], vt2)); 278 tmp = -eval[i]; 279 PetscCall(VecAXPY(vt1,tmp,vt2)); 280 PetscCall(VecNorm(vt1, NORM_INFINITY, &norm)); 281 norm = PetscAbs(norm); 282 if (norm > norm_max) norm_max = norm; 283 /* sniff, and bark if necessary */ 284 if (norm > tols[0]) { 285 PetscCall(PetscPrintf(PETSC_COMM_WORLD," residual violation: %" PetscInt_FMT ", resi: %g\n",i, norm)); 286 } 287 } 288 PetscCall(PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max)); 289 break; 290 default: 291 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%" PetscInt_FMT " is not supported \n",cklvl)); 292 } 293 PetscCall(VecDestroy(&vt2)); 294 PetscCall(VecDestroy(&vt1)); 295 PetscFunctionReturn(0); 296 } 297 298 /*TEST 299 300 build: 301 requires: complex 302 303 test: 304 305 test: 306 suffix: 2 307 args: -test_zheevx 308 309 test: 310 suffix: 3 311 args: -test_zhegv 312 313 test: 314 suffix: 4 315 args: -test_zhegvx 316 317 TEST*/ 318