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