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 PetscErrorCode ierr; 15 PetscBool isSymmetric; 16 PetscScalar *arrayA,*arrayB,*evecs_array=NULL,*work; 17 PetscReal *evals,*rwork; 18 PetscMPIInt size; 19 PetscInt m,i,j,cklvl=2; 20 PetscReal vl,vu,abstol=1.e-8; 21 PetscBLASInt nn,nevs,il,iu,*iwork,*ifail,lwork,lierr,bn,one=1; 22 PetscReal tols[2]; 23 PetscScalar v,sigma2; 24 PetscRandom rctx; 25 PetscReal h2,sigma1 = 100.0; 26 PetscInt dim,Ii,J,n = 6,use_random; 27 28 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 29 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 30 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 31 32 ierr = PetscOptionsHasName(NULL,NULL, "-test_zheevx", &flg);CHKERRQ(ierr); 33 if (flg) { 34 TestZHEEV = PETSC_FALSE; 35 TestZHEEVX = PETSC_TRUE; 36 } 37 ierr = PetscOptionsHasName(NULL,NULL, "-test_zhegv", &flg);CHKERRQ(ierr); 38 if (flg) { 39 TestZHEEV = PETSC_FALSE; 40 TestZHEGV = PETSC_TRUE; 41 } 42 ierr = PetscOptionsHasName(NULL,NULL, "-test_zhegvx", &flg);CHKERRQ(ierr); 43 if (flg) { 44 TestZHEEV = PETSC_FALSE; 45 TestZHEGVX = PETSC_TRUE; 46 } 47 48 ierr = PetscOptionsGetReal(NULL,NULL,"-sigma1",&sigma1,NULL);CHKERRQ(ierr); 49 ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr); 50 dim = n*n; 51 52 ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr); 53 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);CHKERRQ(ierr); 54 ierr = MatSetType(A,MATSEQDENSE);CHKERRQ(ierr); 55 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 56 ierr = MatSetUp(A);CHKERRQ(ierr); 57 58 ierr = PetscOptionsHasName(NULL,NULL,"-norandom",&flg);CHKERRQ(ierr); 59 if (flg) use_random = 0; 60 else use_random = 1; 61 if (use_random) { 62 ierr = PetscRandomCreate(PETSC_COMM_SELF,&rctx);CHKERRQ(ierr); 63 ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr); 64 ierr = PetscRandomSetInterval(rctx,0.0,PETSC_i);CHKERRQ(ierr); 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; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr); 73 } 74 if (i<n-1) { 75 J = Ii+n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr); 76 } 77 if (j>0) { 78 J = Ii-1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr); 79 } 80 if (j<n-1) { 81 J = Ii+1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr); 82 } 83 if (use_random) {ierr = PetscRandomGetValue(rctx,&sigma2);CHKERRQ(ierr);} 84 v = 4.0 - sigma1*h2; 85 ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr); 86 } 87 /* make A complex Hermitian */ 88 v = sigma2*h2; 89 Ii = 0; J = 1; 90 ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr); 91 v = -sigma2*h2; 92 ierr = MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr); 93 if (use_random) {ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);} 94 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 95 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 96 m = n = dim; 97 98 /* Check whether A is symmetric */ 99 ierr = PetscOptionsHasName(NULL,NULL, "-check_symmetry", &flg);CHKERRQ(ierr); 100 if (flg) { 101 Mat Trans; 102 ierr = MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);CHKERRQ(ierr); 103 ierr = MatEqual(A, Trans, &isSymmetric);CHKERRQ(ierr); 104 if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric"); 105 ierr = MatDestroy(&Trans);CHKERRQ(ierr); 106 } 107 108 /* Convert aij matrix to MatSeqDense for LAPACK */ 109 ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 110 if (flg) { 111 ierr = MatDuplicate(A,MAT_COPY_VALUES,&A_dense);CHKERRQ(ierr); 112 } else { 113 ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);CHKERRQ(ierr); 114 } 115 116 ierr = MatCreate(PETSC_COMM_SELF,&B);CHKERRQ(ierr); 117 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,dim,dim);CHKERRQ(ierr); 118 ierr = MatSetType(B,MATSEQDENSE);CHKERRQ(ierr); 119 ierr = MatSetFromOptions(B);CHKERRQ(ierr); 120 ierr = MatSetUp(B);CHKERRQ(ierr); 121 v = 1.0; 122 for (Ii=0; Ii<dim; Ii++) { 123 ierr = MatSetValues(B,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr); 124 } 125 126 /* Solve standard eigenvalue problem: A*x = lambda*x */ 127 /*===================================================*/ 128 ierr = PetscBLASIntCast(2*n,&lwork);CHKERRQ(ierr); 129 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 130 ierr = PetscMalloc1(n,&evals);CHKERRQ(ierr); 131 ierr = PetscMalloc1(lwork,&work);CHKERRQ(ierr); 132 ierr = MatDenseGetArray(A_dense,&arrayA);CHKERRQ(ierr); 133 134 if (TestZHEEV) { /* test zheev() */ 135 ierr = PetscPrintf(PETSC_COMM_WORLD," LAPACKsyev: compute all %D eigensolutions...\n",m);CHKERRQ(ierr); 136 ierr = PetscMalloc1(3*n-2,&rwork);CHKERRQ(ierr); 137 LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,rwork,&lierr); 138 ierr = PetscFree(rwork);CHKERRQ(ierr); 139 140 evecs_array = arrayA; 141 nevs = m; 142 il =1; iu=m; 143 } 144 if (TestZHEEVX) { 145 il = 1; 146 ierr = PetscBLASIntCast((0.2*m),&iu);CHKERRQ(ierr); 147 ierr = PetscPrintf(PETSC_COMM_WORLD," LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu);CHKERRQ(ierr); 148 ierr = PetscMalloc1(m*n+1,&evecs_array);CHKERRQ(ierr); 149 ierr = PetscMalloc1(7*n+1,&rwork);CHKERRQ(ierr); 150 ierr = PetscMalloc1(5*n+1,&iwork);CHKERRQ(ierr); 151 ierr = PetscMalloc1(n+1,&ifail);CHKERRQ(ierr); 152 153 /* in the case "I", vl and vu are not referenced */ 154 vl = 0.0; vu = 8.0; 155 ierr = PetscBLASIntCast(n,&nn);CHKERRQ(ierr); 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 ierr = PetscFree(iwork);CHKERRQ(ierr); 158 ierr = PetscFree(ifail);CHKERRQ(ierr); 159 ierr = PetscFree(rwork);CHKERRQ(ierr); 160 } 161 if (TestZHEGV) { 162 ierr = PetscPrintf(PETSC_COMM_WORLD," LAPACKsygv: compute all %D eigensolutions...\n",m);CHKERRQ(ierr); 163 ierr = PetscMalloc1(3*n+1,&rwork);CHKERRQ(ierr); 164 ierr = MatDenseGetArray(B,&arrayB);CHKERRQ(ierr); 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 ierr = MatDenseRestoreArray(B,&arrayB);CHKERRQ(ierr); 170 ierr = PetscFree(rwork);CHKERRQ(ierr); 171 } 172 if (TestZHEGVX) { 173 il = 1; 174 ierr = PetscBLASIntCast((0.2*m),&iu);CHKERRQ(ierr); 175 ierr = PetscPrintf(PETSC_COMM_WORLD," LAPACKsygv: compute %d to %d-th eigensolutions...\n",il,iu);CHKERRQ(ierr); 176 ierr = PetscMalloc1(m*n+1,&evecs_array);CHKERRQ(ierr); 177 ierr = PetscMalloc1(6*n+1,&iwork);CHKERRQ(ierr); 178 ifail = iwork + 5*n; 179 ierr = PetscMalloc1(7*n+1,&rwork);CHKERRQ(ierr); 180 ierr = MatDenseGetArray(B,&arrayB);CHKERRQ(ierr); 181 vl = 0.0; vu = 8.0; 182 ierr = PetscBLASIntCast(n,&nn);CHKERRQ(ierr); 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 ierr = MatDenseRestoreArray(B,&arrayB);CHKERRQ(ierr); 185 ierr = PetscFree(iwork);CHKERRQ(ierr); 186 ierr = PetscFree(rwork);CHKERRQ(ierr); 187 } 188 ierr = MatDenseRestoreArray(A_dense,&arrayA);CHKERRQ(ierr); 189 if (nevs <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%d, no eigensolution has found", nevs); 190 191 /* View evals */ 192 ierr = PetscOptionsHasName(NULL,NULL, "-eig_view", &flg);CHKERRQ(ierr); 193 if (flg) { 194 ierr = PetscPrintf(PETSC_COMM_WORLD," %d evals: \n",nevs);CHKERRQ(ierr); 195 for (i=0; i<nevs; i++) {ierr = PetscPrintf(PETSC_COMM_WORLD,"%D %g\n",i+il,(double)evals[i]);CHKERRQ(ierr);} 196 } 197 198 /* Check residuals and orthogonality */ 199 ierr = PetscMalloc1(nevs+1,&evecs);CHKERRQ(ierr); 200 for (i=0; i<nevs; i++) { 201 ierr = VecCreate(PETSC_COMM_SELF,&evecs[i]);CHKERRQ(ierr); 202 ierr = VecSetSizes(evecs[i],PETSC_DECIDE,n);CHKERRQ(ierr); 203 ierr = VecSetFromOptions(evecs[i]);CHKERRQ(ierr); 204 ierr = VecPlaceArray(evecs[i],evecs_array+i*n);CHKERRQ(ierr); 205 } 206 207 tols[0] = PETSC_SQRT_MACHINE_EPSILON; tols[1] = PETSC_SQRT_MACHINE_EPSILON; 208 ierr = CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);CHKERRQ(ierr); 209 for (i=0; i<nevs; i++) { ierr = VecDestroy(&evecs[i]);CHKERRQ(ierr);} 210 ierr = PetscFree(evecs);CHKERRQ(ierr); 211 212 /* Free work space. */ 213 if (TestZHEEVX || TestZHEGVX) { 214 ierr = PetscFree(evecs_array);CHKERRQ(ierr); 215 } 216 ierr = PetscFree(evals);CHKERRQ(ierr); 217 ierr = PetscFree(work);CHKERRQ(ierr); 218 ierr = MatDestroy(&A_dense);CHKERRQ(ierr); 219 ierr = MatDestroy(&A);CHKERRQ(ierr); 220 ierr = MatDestroy(&B);CHKERRQ(ierr); 221 ierr = PetscFinalize(); 222 return ierr; 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 ierr,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 ierr = VecDuplicate(evec[0],&vt1);CHKERRQ(ierr); 250 ierr = VecDuplicate(evec[0],&vt2);CHKERRQ(ierr); 251 252 switch (cklvl) { 253 case 2: 254 dot_max = 0.0; 255 for (i = il; i<iu; i++) { 256 ierr = VecCopy(evec[i], vt1);CHKERRQ(ierr); 257 for (j=il; j<iu; j++) { 258 ierr = VecDot(evec[j],vt1,&dot);CHKERRQ(ierr); 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 ierr = VecNorm(evec[i],NORM_INFINITY,&norm);CHKERRQ(ierr); 267 ierr = PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)rdot,(double)norm);CHKERRQ(ierr); 268 } 269 } 270 } 271 ierr = PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);CHKERRQ(ierr); 272 273 case 1: 274 norm_max = 0.0; 275 for (i = il; i< iu; i++) { 276 ierr = MatMult(A, evec[i], vt1);CHKERRQ(ierr); 277 ierr = VecCopy(evec[i], vt2);CHKERRQ(ierr); 278 tmp = -eval[i]; 279 ierr = VecAXPY(vt1,tmp,vt2);CHKERRQ(ierr); 280 ierr = VecNorm(vt1, NORM_INFINITY, &norm);CHKERRQ(ierr); 281 norm = PetscAbs(norm); 282 if (norm > norm_max) norm_max = norm; 283 /* sniff, and bark if necessary */ 284 if (norm > tols[0]) { 285 ierr = PetscPrintf(PETSC_COMM_WORLD," residual violation: %d, resi: %g\n",i, norm);CHKERRQ(ierr); 286 } 287 } 288 ierr = PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max);CHKERRQ(ierr); 289 break; 290 default: 291 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);CHKERRQ(ierr); 292 } 293 ierr = VecDestroy(&vt2);CHKERRQ(ierr); 294 ierr = VecDestroy(&vt1);CHKERRQ(ierr); 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