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


 10: int main(int argc,char **args)
 11: {
 12: #if !defined(PETSC_USE_COMPLEX)
 14:   PetscReal      *work,tols[2];
 15:   PetscInt       i,j,n,il=1,iu=5,*iblock,*isplit,*iwork,nevs,*ifail,cklvl=2;
 16:   PetscMPIInt    size;
 17:   PetscTruth     flg;
 18:   Vec            *evecs;
 19:   PetscScalar    *evecs_array,*D,*E,*evals;
 20:   Mat            T;
 21: #if !defined(PETSC_MISSING_LAPACK_DSTEBZ)
 22:   PetscReal      vl=0.0,vu=4.0,tol=1.e-10;
 23:   PetscInt       nsplit,info;
 24: #endif
 25: #endif
 26: 
 27:   PetscInitialize(&argc,&args,(char *)0,help);
 28: #if defined(PETSC_USE_COMPLEX)
 29:   SETERRQ(1,"This example does not work with complex numbers");
 30: #else
 31:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 32:   if (size != 1) SETERRQ(PETSC_ERR_SUP,"This is a uniprocessor example only!");

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

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

 50:   /* Solve eigenvalue problem: A*evec = eval*evec */
 51: #if defined(PETSC_MISSING_LAPACK_STEBZ)
 52:   SETERRQ(PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
 53: #else
 54:   printf(" LAPACKstebz_: compute %d eigenvalues...\n",nevs);
 55:   LAPACKstebz_("I","E",&n,&vl,&vu,&il,&iu,&tol,(PetscReal*)D,(PetscReal*)E,&nevs,&nsplit,(PetscReal*)evals,iblock,isplit,work,iwork,&info);
 56:   if (info) SETERRQ1(PETSC_ERR_USER,"LAPACKstebz_ fails. info %d",info);
 57: #endif

 59:   printf(" LAPACKstein_: compute %d found eigenvectors...\n",nevs);
 60:   PetscMalloc(n*nevs*sizeof(PetscScalar),&evecs_array);
 61:   PetscMalloc(nevs*sizeof(PetscInt),&ifail);
 62: #if defined(PETSC_MISSING_LAPACK_STEIN)
 63:   SETERRQ(PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
 64: #else
 65:   LAPACKstein_(&n,(PetscReal*)D,(PetscReal*)E,&nevs,(PetscReal*)evals,iblock,isplit,evecs_array,&n,work,iwork,ifail,&info);
 66:   if (info) SETERRQ1(PETSC_ERR_USER,"LAPACKstein_ fails. info %d",info);
 67: #endif
 68:   /* View evals */
 69:   PetscOptionsHasName(PETSC_NULL, "-eig_view", &flg);
 70:   if (flg){
 71:     printf(" %d evals: \n",nevs);
 72:     for (i=0; i<nevs; i++) printf("%d  %G\n",i,evals[i]);
 73:   }

 75:   /* Check residuals and orthogonality */
 76:   MatCreate(PETSC_COMM_SELF,&T);
 77:   MatSetSizes(T,PETSC_DECIDE,PETSC_DECIDE,n,n);
 78:   MatSetType(T,MATSBAIJ);
 79:   MatSetFromOptions(T);
 80:   for (i=0; i<n; i++){
 81:     MatSetValues(T,1,&i,1,&i,&D[i],INSERT_VALUES);
 82:     if (i != n-1){
 83:       j = i+1;
 84:       MatSetValues(T,1,&i,1,&j,&E[i],INSERT_VALUES);
 85:     }
 86:   }
 87:   MatAssemblyBegin(T,MAT_FINAL_ASSEMBLY);
 88:   MatAssemblyEnd(T,MAT_FINAL_ASSEMBLY);

 90:   PetscMalloc((nevs+1)*sizeof(Vec),&evecs);
 91:   for (i=0; i<nevs; i++){
 92:     VecCreate(PETSC_COMM_SELF,&evecs[i]);
 93:     VecSetSizes(evecs[i],PETSC_DECIDE,n);
 94:     VecSetFromOptions(evecs[i]);
 95:     VecPlaceArray(evecs[i],evecs_array+i*n);
 96:   }
 97: 
 98:   tols[0] = 1.e-8;  tols[1] = 1.e-8;
 99:   CkEigenSolutions(cklvl,T,il-1,iu-1,evals,evecs,tols);

101:   for (i=0; i<nevs; i++){
102:     VecResetArray(evecs[i]);
103:   }

105:   /* free space */

107:   MatDestroy(T);

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

147:   nev = iu - il;
148:   if (nev <= 0) return(0);

150:   VecDuplicate(evec[0],&vt1);
151:   VecDuplicate(evec[0],&vt2);

153:   switch (cklvl){
154:   case 2:
155:     dot_max = 0.0;
156:     for (i = il; i<iu; i++){
157:       VecCopy(evec[i], vt1);
158:       for (j=il; j<iu; j++){
159:         VecDot(evec[j],vt1,&dot);
160:         if (j == i){
161:           dot = PetscAbsScalar(dot - 1.0);
162:         } else {
163:           dot = PetscAbsScalar(dot);
164:         }
165:         if (PetscAbsScalar(dot) > dot_max) dot_max = PetscAbsScalar(dot);
166: #ifdef DEBUG_CkEigenSolutions
167:         if (dot > tols[1] ) {
168:           VecNorm(evec[i],NORM_INFINITY,&norm);
169:           PetscPrintf(PETSC_COMM_SELF,"|delta(%d,%d)|: %G, norm: %G\n",i,j,dot,norm);
170:         }
171: #endif
172:       }
173:     }
174:     PetscPrintf(PETSC_COMM_SELF,"    max|(x_j^T*x_i) - delta_ji|: %G\n",dot_max);

176:   case 1:
177:     norm_max = 0.0;
178:     for (i = il; i< iu; i++){
179:       MatMult(A, evec[i], vt1);
180:       VecCopy(evec[i], vt2);
181:       tmp  = -eval[i];
182:       VecAXPY(vt1,tmp,vt2);
183:       VecNorm(vt1, NORM_INFINITY, &norm);
184:       norm = PetscAbsScalar(norm);
185:       if (norm > norm_max) norm_max = norm;
186: #ifdef DEBUG_CkEigenSolutions
187:       /* sniff, and bark if necessary */
188:       if (norm > tols[0]){
189:         printf( "  residual violation: %d, resi: %g\n",i, norm);
190:       }
191: #endif
192:     }
193:     PetscPrintf(PETSC_COMM_SELF,"    max_resi:                    %G\n", norm_max);
194:    break;
195:   default:
196:     PetscPrintf(PETSC_COMM_SELF,"Error: cklvl=%d is not supported \n",cklvl);
197:   }

199:   VecDestroy(vt2);
200:   VecDestroy(vt1);
201:   return(0);
202: }