Actual source code: ex118.c

  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*);

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

 25:   PetscInitialize(&argc,&args,(char*)0,help);
 26:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 29:   n      = 100;
 30:   nevs   = iu - il;
 31:   PetscMalloc1(3*n+1,&D);
 32:   E      = D + n;
 33:   evals  = E + n;
 34:   PetscMalloc1(5*n+1,&work);
 35:   PetscMalloc1(3*n+1,&iwork);
 36:   PetscMalloc1(3*n+1,&iblock);
 37:   isplit = iblock + n;

 39:   /* Set symmetric tridiagonal matrix */
 40:   for (i=0; i<n; i++) {
 41:     D[i] = 2.0;
 42:     E[i] = 1.0;
 43:   }

 45:   /* Solve eigenvalue problem: A*evec = eval*evec */
 46:   PetscPrintf(PETSC_COMM_SELF," LAPACKstebz_: compute %d eigenvalues...\n",nevs);
 47:   LAPACKstebz_("I","E",&n,&vl,&vu,&il,&iu,&tol,(PetscReal*)D,(PetscReal*)E,&nevs,&nsplit,(PetscReal*)evals,iblock,isplit,work,iwork,&info);

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

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

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

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

 89:   for (i=0; i<nevs; i++) {
 90:     VecResetArray(evecs[i]);
 91:   }

 93:   /* free space */

 95:   MatDestroy(&T);

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

132:   nev = iu - il;
133:   if (nev <= 0) return 0;

135:   VecDuplicate(evec[0],&vt1);
136:   VecDuplicate(evec[0],&vt2);

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

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

183:   VecDestroy(&vt2);
184:   VecDestroy(&vt1);
185:   return 0;
186: }