1 static char help[] = "Test LAPACK routine DSYEV() or DSYEVX(). \n\ 2 Reads PETSc matrix A \n\ 3 then computes selected eigenvalues, and optionally, eigenvectors of \n\ 4 a real generalized symmetric-definite eigenproblem \n\ 5 A*x = lambda*x \n\ 6 Input parameters include\n\ 7 -f <input_file> : file to load\n\ 8 e.g. ./ex116 -f $DATAFILESPATH/matrices/small \n\n"; 9 10 #include <petscmat.h> 11 #include <petscblaslapack.h> 12 13 extern PetscErrorCode CkEigenSolutions(PetscInt,Mat,PetscInt,PetscInt,PetscReal*,Vec*,PetscReal*); 14 15 int main(int argc,char **args) 16 { 17 Mat A,A_dense; 18 Vec *evecs; 19 PetscViewer fd; /* viewer */ 20 char file[1][PETSC_MAX_PATH_LEN]; /* input file name */ 21 PetscBool flg,TestSYEVX=PETSC_TRUE; 22 PetscErrorCode ierr; 23 PetscBool isSymmetric; 24 PetscScalar *arrayA,*evecs_array,*work,*evals; 25 PetscMPIInt size; 26 PetscInt m,n,i,cklvl=2; 27 PetscBLASInt nevs,il,iu,in; 28 PetscReal vl,vu,abstol=1.e-8; 29 PetscBLASInt *iwork,*ifail,lwork,lierr,bn; 30 PetscReal tols[2]; 31 32 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 33 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 34 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 35 36 ierr = PetscOptionsHasName(NULL,NULL, "-test_syev", &flg);CHKERRQ(ierr); 37 if (flg) { 38 TestSYEVX = PETSC_FALSE; 39 } 40 41 /* Determine files from which we read the two matrices */ 42 ierr = PetscOptionsGetString(NULL,NULL,"-f",file[0],sizeof(file[0]),&flg);CHKERRQ(ierr); 43 44 /* Load matrix A */ 45 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&fd);CHKERRQ(ierr); 46 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 47 ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr); 48 ierr = MatLoad(A,fd);CHKERRQ(ierr); 49 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 50 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 51 52 /* Check whether A is symmetric */ 53 ierr = PetscOptionsHasName(NULL,NULL, "-check_symmetry", &flg);CHKERRQ(ierr); 54 if (flg) { 55 Mat Trans; 56 ierr = MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);CHKERRQ(ierr); 57 ierr = MatEqual(A, Trans, &isSymmetric);CHKERRQ(ierr); 58 if (!isSymmetric) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A must be symmetric"); 59 ierr = MatDestroy(&Trans);CHKERRQ(ierr); 60 } 61 62 /* Solve eigenvalue problem: A_dense*x = lambda*B*x */ 63 /*==================================================*/ 64 /* Convert aij matrix to MatSeqDense for LAPACK */ 65 ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);CHKERRQ(ierr); 66 67 ierr = PetscBLASIntCast(8*n,&lwork);CHKERRQ(ierr); 68 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 69 ierr = PetscMalloc1(n,&evals);CHKERRQ(ierr); 70 ierr = PetscMalloc1(lwork,&work);CHKERRQ(ierr); 71 ierr = MatDenseGetArray(A_dense,&arrayA);CHKERRQ(ierr); 72 73 if (!TestSYEVX) { /* test syev() */ 74 ierr = PetscPrintf(PETSC_COMM_SELF," LAPACKsyev: compute all %D eigensolutions...\n",m);CHKERRQ(ierr); 75 LAPACKsyev_("V","U",&bn,arrayA,&bn,evals,work,&lwork,&lierr); 76 evecs_array = arrayA; 77 ierr = PetscBLASIntCast(m,&nevs);CHKERRQ(ierr); 78 il = 1; 79 ierr = PetscBLASIntCast(m,&iu);CHKERRQ(ierr); 80 } else { /* test syevx() */ 81 il = 1; 82 ierr = PetscBLASIntCast(0.2*m,&iu);CHKERRQ(ierr); 83 ierr = PetscBLASIntCast(n,&in);CHKERRQ(ierr); 84 ierr = PetscPrintf(PETSC_COMM_SELF," LAPACKsyevx: compute %d to %d-th eigensolutions...\n",il,iu);CHKERRQ(ierr); 85 ierr = PetscMalloc1(m*n+1,&evecs_array);CHKERRQ(ierr); 86 ierr = PetscMalloc1(6*n+1,&iwork);CHKERRQ(ierr); 87 ifail = iwork + 5*n; 88 89 /* in the case "I", vl and vu are not referenced */ 90 vl = 0.0; vu = 8.0; 91 LAPACKsyevx_("V","I","U",&bn,arrayA,&bn,&vl,&vu,&il,&iu,&abstol,&nevs,evals,evecs_array,&in,work,&lwork,iwork,ifail,&lierr); 92 ierr = PetscFree(iwork);CHKERRQ(ierr); 93 } 94 ierr = MatDenseRestoreArray(A_dense,&arrayA);CHKERRQ(ierr); 95 if (nevs <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED, "nev=%D, no eigensolution has found", nevs); 96 97 /* View eigenvalues */ 98 ierr = PetscOptionsHasName(NULL,NULL, "-eig_view", &flg);CHKERRQ(ierr); 99 if (flg) { 100 ierr = PetscPrintf(PETSC_COMM_SELF," %D evals: \n",nevs);CHKERRQ(ierr); 101 for (i=0; i<nevs; i++) {ierr = PetscPrintf(PETSC_COMM_SELF,"%D %g\n",i+il,(double)evals[i]);CHKERRQ(ierr);} 102 } 103 104 /* Check residuals and orthogonality */ 105 ierr = PetscMalloc1(nevs+1,&evecs);CHKERRQ(ierr); 106 for (i=0; i<nevs; i++) { 107 ierr = VecCreate(PETSC_COMM_SELF,&evecs[i]);CHKERRQ(ierr); 108 ierr = VecSetSizes(evecs[i],PETSC_DECIDE,n);CHKERRQ(ierr); 109 ierr = VecSetFromOptions(evecs[i]);CHKERRQ(ierr); 110 ierr = VecPlaceArray(evecs[i],evecs_array+i*n);CHKERRQ(ierr); 111 } 112 113 tols[0] = tols[1] = PETSC_SQRT_MACHINE_EPSILON; 114 ierr = CkEigenSolutions(cklvl,A,il-1,iu-1,evals,evecs,tols);CHKERRQ(ierr); 115 116 /* Free work space. */ 117 for (i=0; i<nevs; i++) { ierr = VecDestroy(&evecs[i]);CHKERRQ(ierr);} 118 ierr = PetscFree(evecs);CHKERRQ(ierr); 119 ierr = MatDestroy(&A_dense);CHKERRQ(ierr); 120 ierr = PetscFree(work);CHKERRQ(ierr); 121 if (TestSYEVX) {ierr = PetscFree(evecs_array);CHKERRQ(ierr);} 122 123 /* Compute SVD: A_dense = U*SIGMA*transpose(V), 124 JOBU=JOBV='S': the first min(m,n) columns of U and V are returned in the arrayU and arrayV; */ 125 /*==============================================================================================*/ 126 { 127 /* Convert aij matrix to MatSeqDense for LAPACK */ 128 PetscScalar *arrayU,*arrayVT,*arrayErr,alpha=1.0,beta=-1.0; 129 Mat Err; 130 PetscBLASInt minMN,maxMN,im,in; 131 PetscInt j; 132 PetscReal norm; 133 134 ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&A_dense);CHKERRQ(ierr); 135 136 minMN = PetscMin(m,n); 137 maxMN = PetscMax(m,n); 138 lwork = 5*minMN + maxMN; 139 ierr = PetscMalloc4(m*minMN,&arrayU,m*minMN,&arrayVT,m*minMN,&arrayErr,lwork,&work);CHKERRQ(ierr); 140 141 /* Create matrix Err for checking error */ 142 ierr = MatCreate(PETSC_COMM_WORLD,&Err);CHKERRQ(ierr); 143 ierr = MatSetSizes(Err,PETSC_DECIDE,PETSC_DECIDE,m,minMN);CHKERRQ(ierr); 144 ierr = MatSetType(Err,MATSEQDENSE);CHKERRQ(ierr); 145 ierr = MatSeqDenseSetPreallocation(Err,(PetscScalar*)arrayErr);CHKERRQ(ierr); 146 147 /* Save A to arrayErr for checking accuracy later. arrayA will be destroyed by LAPACKgesvd_() */ 148 ierr = MatDenseGetArray(A_dense,&arrayA);CHKERRQ(ierr); 149 ierr = PetscArraycpy(arrayErr,arrayA,m*minMN);CHKERRQ(ierr); 150 151 ierr = PetscBLASIntCast(m,&im);CHKERRQ(ierr); 152 ierr = PetscBLASIntCast(n,&in);CHKERRQ(ierr); 153 /* Compute A = U*SIGMA*VT */ 154 LAPACKgesvd_("S","S",&im,&in,arrayA,&im,evals,arrayU,&minMN,arrayVT,&minMN,work,&lwork,&lierr); 155 ierr = MatDenseRestoreArray(A_dense,&arrayA);CHKERRQ(ierr); 156 if (!lierr) { 157 ierr = PetscPrintf(PETSC_COMM_SELF," 1st 10 of %d singular values: \n",minMN);CHKERRQ(ierr); 158 for (i=0; i<10; i++) {ierr = PetscPrintf(PETSC_COMM_SELF,"%D %g\n",i,(double)evals[i]);CHKERRQ(ierr);} 159 } else { 160 ierr = PetscPrintf(PETSC_COMM_SELF,"LAPACKgesvd_ fails!");CHKERRQ(ierr); 161 } 162 163 /* Check Err = (U*Sigma*V^T - A) using BLASgemm() */ 164 /* U = U*Sigma */ 165 for (j=0; j<minMN; j++) { /* U[:,j] = sigma[j]*U[:,j] */ 166 for (i=0; i<m; i++) arrayU[j*m+i] *= evals[j]; 167 } 168 /* Err = U*VT - A = alpha*U*VT + beta*Err */ 169 BLASgemm_("N","N",&im,&minMN,&minMN,&alpha,arrayU,&im,arrayVT,&minMN,&beta,arrayErr,&im); 170 ierr = MatNorm(Err,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 171 ierr = PetscPrintf(PETSC_COMM_SELF," || U*Sigma*VT - A || = %g\n",(double)norm);CHKERRQ(ierr); 172 173 ierr = PetscFree4(arrayU,arrayVT,arrayErr,work);CHKERRQ(ierr); 174 ierr = PetscFree(evals);CHKERRQ(ierr); 175 ierr = MatDestroy(&A_dense);CHKERRQ(ierr); 176 ierr = MatDestroy(&Err);CHKERRQ(ierr); 177 } 178 179 ierr = MatDestroy(&A);CHKERRQ(ierr); 180 ierr = PetscFinalize(); 181 return ierr; 182 } 183 /*------------------------------------------------ 184 Check the accuracy of the eigen solution 185 ----------------------------------------------- */ 186 /* 187 input: 188 cklvl - check level: 189 1: check residual 190 2: 1 and check B-orthogonality locally 191 A - matrix 192 il,iu - lower and upper index bound of eigenvalues 193 eval, evec - eigenvalues and eigenvectors stored in this process 194 tols[0] - reporting tol_res: || A * evec[i] - eval[i]*evec[i] || 195 tols[1] - reporting tol_orth: evec[i]^T*evec[j] - delta_ij 196 */ 197 PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscReal *eval,Vec *evec,PetscReal *tols) 198 { 199 PetscInt ierr,i,j,nev; 200 Vec vt1,vt2; /* tmp vectors */ 201 PetscReal norm,tmp,dot,norm_max,dot_max; 202 203 PetscFunctionBegin; 204 nev = iu - il; 205 if (nev <= 0) PetscFunctionReturn(0); 206 207 /*ierr = VecView(evec[0],PETSC_VIEWER_STDOUT_WORLD);*/ 208 ierr = VecDuplicate(evec[0],&vt1);CHKERRQ(ierr); 209 ierr = VecDuplicate(evec[0],&vt2);CHKERRQ(ierr); 210 211 switch (cklvl) { 212 case 2: 213 dot_max = 0.0; 214 for (i = il; i<iu; i++) { 215 ierr = VecCopy(evec[i], vt1);CHKERRQ(ierr); 216 for (j=il; j<iu; j++) { 217 ierr = VecDot(evec[j],vt1,&dot);CHKERRQ(ierr); 218 if (j == i) { 219 dot = PetscAbsScalar(dot - 1); 220 } else { 221 dot = PetscAbsScalar(dot); 222 } 223 if (dot > dot_max) dot_max = dot; 224 if (dot > tols[1]) { 225 ierr = VecNorm(evec[i],NORM_INFINITY,&norm);CHKERRQ(ierr); 226 ierr = PetscPrintf(PETSC_COMM_SELF,"|delta(%D,%D)|: %g, norm: %g\n",i,j,(double)dot,(double)norm);CHKERRQ(ierr); 227 } 228 } 229 } 230 ierr = PetscPrintf(PETSC_COMM_SELF," max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);CHKERRQ(ierr); 231 232 case 1: 233 norm_max = 0.0; 234 for (i = il; i< iu; i++) { 235 ierr = MatMult(A, evec[i], vt1);CHKERRQ(ierr); 236 ierr = VecCopy(evec[i], vt2);CHKERRQ(ierr); 237 tmp = -eval[i]; 238 ierr = VecAXPY(vt1,tmp,vt2);CHKERRQ(ierr); 239 ierr = VecNorm(vt1, NORM_INFINITY, &norm);CHKERRQ(ierr); 240 norm = PetscAbsScalar(norm); 241 if (norm > norm_max) norm_max = norm; 242 /* sniff, and bark if necessary */ 243 if (norm > tols[0]) { 244 ierr = PetscPrintf(PETSC_COMM_SELF," residual violation: %D, resi: %g\n",i, norm);CHKERRQ(ierr); 245 } 246 } 247 ierr = PetscPrintf(PETSC_COMM_SELF," max_resi: %g\n", (double)norm_max);CHKERRQ(ierr); 248 break; 249 default: 250 ierr = PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%D is not supported \n",cklvl);CHKERRQ(ierr); 251 } 252 ierr = VecDestroy(&vt2);CHKERRQ(ierr); 253 ierr = VecDestroy(&vt1);CHKERRQ(ierr); 254 PetscFunctionReturn(0); 255 } 256 257 /*TEST 258 259 build: 260 requires: !complex 261 262 test: 263 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 264 args: -f ${DATAFILESPATH}/matrices/small 265 output_file: output/ex116_1.out 266 267 test: 268 suffix: 2 269 requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES) 270 args: -f ${DATAFILESPATH}/matrices/small -test_syev -check_symmetry 271 272 TEST*/ 273