Actual source code: ex118.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: static char help[] = "Test LAPACK routine DSTEBZ() and DTEIN().  \n\n";

  3: #include <petscmat.h>
  4: #include <petscblaslapack.h>

  6: extern PetscErrorCode CkEigenSolutions(PetscInt,Mat,PetscInt,PetscInt,PetscScalar*,Vec*,PetscReal*);

 10: int main(int argc,char **args)
 11: {
 12: #if defined(PETSC_USE_COMPLEX) || defined(PETSC_MISSING_LAPACK_DSTEBZ) || defined(PETSC_MISSING_LAPACK_STEIN)
 13:   PetscInitialize(&argc,&args,(char*)0,help);
 14:   SETERRQ(PETSC_COMM_WORLD,1,"This example requires LAPACK routines dstebz and stien and real numbers");
 15: #else
 17:   PetscReal      *work,tols[2];
 18:   PetscInt       i,j;
 19:   PetscBLASInt   n,il=1,iu=5,*iblock,*isplit,*iwork,nevs,*ifail,cklvl=2;
 20:   PetscMPIInt    size;
 21:   PetscBool      flg;
 22:   Vec            *evecs;
 23:   PetscScalar    *evecs_array,*D,*E,*evals;
 24:   Mat            T;
 25:   PetscReal      vl=0.0,vu=4.0,tol= 1000*PETSC_MACHINE_EPSILON;
 26:   PetscBLASInt   nsplit,info;


 29:   PetscInitialize(&argc,&args,(char*)0,help);
 30:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 31:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");

 33:   n      = 100;
 34:   nevs   = iu - il;
 35:   PetscMalloc1(3*n+1,&D);
 36:   E      = D + n;
 37:   evals  = E + n;
 38:   PetscMalloc1(5*n+1,&work);
 39:   PetscMalloc1(3*n+1,&iwork);
 40:   PetscMalloc1(3*n+1,&iblock);
 41:   isplit = iblock + n;

 43:   /* Set symmetric tridiagonal matrix */
 44:   for (i=0; i<n; i++) {
 45:     D[i] = 2.0;
 46:     E[i] = 1.0;
 47:   }

 49:   /* Solve eigenvalue problem: A*evec = eval*evec */
 50:   PetscPrintf(PETSC_COMM_SELF," LAPACKstebz_: compute %d eigenvalues...\n",nevs);
 51:   LAPACKstebz_("I","E",&n,&vl,&vu,&il,&iu,&tol,(PetscReal*)D,(PetscReal*)E,&nevs,&nsplit,(PetscReal*)evals,iblock,isplit,work,iwork,&info);
 52:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstebz_ fails. info %d",info);

 54:   PetscPrintf(PETSC_COMM_SELF," LAPACKstein_: compute %d found eigenvectors...\n",nevs);
 55:   PetscMalloc1(n*nevs,&evecs_array);
 56:   PetscMalloc1(nevs,&ifail);
 57:   LAPACKstein_(&n,(PetscReal*)D,(PetscReal*)E,&nevs,(PetscReal*)evals,iblock,isplit,evecs_array,&n,work,iwork,ifail,&info);
 58:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"LAPACKstein_ fails. info %d",info);
 59:   /* View evals */
 60:   PetscOptionsHasName(NULL,NULL, "-eig_view", &flg);
 61:   if (flg) {
 62:     PetscPrintf(PETSC_COMM_SELF," %d evals: \n",nevs);
 63:     for (i=0; i<nevs; i++) {PetscPrintf(PETSC_COMM_SELF,"%D  %g\n",i,(double)evals[i]);}
 64:   }

 66:   /* Check residuals and orthogonality */
 67:   MatCreate(PETSC_COMM_SELF,&T);
 68:   MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n);
 69:   MatSetType(T,MATSBAIJ);
 70:   MatSetFromOptions(T);
 71:   MatSetUp(T);
 72:   for (i=0; i<n; i++) {
 73:     MatSetValues(T,1,&i,1,&i,&D[i],INSERT_VALUES);
 74:     if (i != n-1) {
 75:       j    = i+1;
 76:       MatSetValues(T,1,&i,1,&j,&E[i],INSERT_VALUES);
 77:     }
 78:   }
 79:   MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);
 80:   MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);

 82:   PetscMalloc1(nevs+1,&evecs);
 83:   for (i=0; i<nevs; i++) {
 84:     VecCreate(PETSC_COMM_SELF,&evecs[i]);
 85:     VecSetSizes(evecs[i],PETSC_DECIDE,n);
 86:     VecSetFromOptions(evecs[i]);
 87:     VecPlaceArray(evecs[i],evecs_array+i*n);
 88:   }

 90:   tols[0] = 1.e-8;  tols[1] = 1.e-8;
 91:   CkEigenSolutions(cklvl,T,il-1,iu-1,evals,evecs,tols);

 93:   for (i=0; i<nevs; i++) {
 94:     VecResetArray(evecs[i]);
 95:   }

 97:   /* free space */

 99:   MatDestroy(&T);

101:   for (i=0; i<nevs; i++) { VecDestroy(&evecs[i]);}
102:   PetscFree(evecs);
103:   PetscFree(D);
104:   PetscFree(work);
105:   PetscFree(iwork);
106:   PetscFree(iblock);
107:   PetscFree(evecs_array);
108:   PetscFree(ifail);
109:   PetscFinalize();
110:   return 0;
111: #endif
112: }
113: /*------------------------------------------------
114:   Check the accuracy of the eigen solution
115:   ----------------------------------------------- */
116: /*
117:   input:
118:      cklvl      - check level:
119:                     1: check residual
120:                     2: 1 and check B-orthogonality locally
121:      A          - matrix
122:      il,iu      - lower and upper index bound of eigenvalues
123:      eval, evec - eigenvalues and eigenvectors stored in this process
124:      tols[0]    - reporting tol_res: || A * evec[i] - eval[i]*evec[i] ||
125:      tols[1]    - reporting tol_orth: evec[i]^T*evec[j] - delta_ij
126: */
127: #undef DEBUG_CkEigenSolutions
130: PetscErrorCode CkEigenSolutions(PetscInt cklvl,Mat A,PetscInt il,PetscInt iu,PetscScalar *eval,Vec *evec,PetscReal *tols)
131: {
132:   PetscInt    ierr,i,j,nev;
133:   Vec         vt1,vt2;  /* tmp vectors */
134:   PetscReal   norm,norm_max;
135:   PetscScalar dot,tmp;
136:   PetscReal   dot_max;

139:   nev = iu - il;
140:   if (nev <= 0) return(0);

142:   VecDuplicate(evec[0],&vt1);
143:   VecDuplicate(evec[0],&vt2);

145:   switch (cklvl) {
146:   case 2:
147:     dot_max = 0.0;
148:     for (i = il; i<iu; i++) {
149:       VecCopy(evec[i], vt1);
150:       for (j=il; j<iu; j++) {
151:         VecDot(evec[j],vt1,&dot);
152:         if (j == i) {
153:           dot = PetscAbsScalar(dot - 1.0);
154:         } else {
155:           dot = PetscAbsScalar(dot);
156:         }
157:         if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot);
158: #if defined(DEBUG_CkEigenSolutions)
159:         if (dot > tols[1]) {
160:           VecNorm(evec[i],NORM_INFINITY,&norm);
161:           PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %g, norm: %d\n",i,j,(double)dot,(double)norm);
162:         }
163: #endif
164:       }
165:     }
166:     PetscPrintf(PETSC_COMM_SELF,"    max|(x_j^T*x_i) - delta_ji|: %g\n",(double)dot_max);

168:   case 1:
169:     norm_max = 0.0;
170:     for (i = il; i< iu; i++) {
171:       MatMult(A, evec[i], vt1);
172:       VecCopy(evec[i], vt2);
173:       tmp  = -eval[i];
174:       VecAXPY(vt1,tmp,vt2);
175:       VecNorm(vt1, NORM_INFINITY, &norm);
176:       norm = PetscAbsReal(norm);
177:       if (norm > norm_max) norm_max = norm;
178: #if defined(DEBUG_CkEigenSolutions)
179:       if (norm > tols[0]) {
180:         PetscPrintf(PETSC_COMM_SELF,"  residual violation: %d, resi: %g\n",i, norm);
181:       }
182: #endif
183:     }
184:     PetscPrintf(PETSC_COMM_SELF,"    max_resi:                    %g\n", (double)norm_max);
185:     break;
186:   default:
187:     PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
188:   }

190:   VecDestroy(&vt2);
191:   VecDestroy(&vt1);
192:   return(0);
193: }