Actual source code: eige.c

  1: /*$Id: eige.c,v 1.35 2001/08/07 03:03:45 balay Exp $*/

 3:  #include src/sles/ksp/kspimpl.h

  5: #undef __FUNCT__  
  7: /*@
  8:     KSPComputeExplicitOperator - Computes the explicit preconditioned operator.  

 10:     Collective on KSP

 12:     Input Parameter:
 13: .   ksp - the Krylov subspace context

 15:     Output Parameter:
 16: .   mat - the explict preconditioned operator

 18:     Notes:
 19:     This computation is done by applying the operators to columns of the 
 20:     identity matrix.

 22:     Currently, this routine uses a dense matrix format when 1 processor
 23:     is used and a sparse format otherwise.  This routine is costly in general,
 24:     and is recommended for use only with relatively small systems.

 26:     Level: advanced
 27:    
 28: .keywords: KSP, compute, explicit, operator

 30: .seealso: KSPComputeEigenvaluesExplicitly()
 31: @*/
 32: int KSPComputeExplicitOperator(KSP ksp,Mat *mat)
 33: {
 34:   Vec      in,out;
 35:   int      ierr,i,M,m,size,*rows,start,end;
 36:   Mat      A;
 37:   MPI_Comm comm;
 38:   PetscScalar   *array,zero = 0.0,one = 1.0;

 43:   comm = ksp->comm;

 45:   MPI_Comm_size(comm,&size);

 47:   VecDuplicate(ksp->vec_sol,&in);
 48:   VecDuplicate(ksp->vec_sol,&out);
 49:   VecGetSize(in,&M);
 50:   VecGetLocalSize(in,&m);
 51:   VecGetOwnershipRange(in,&start,&end);
 52:   PetscMalloc((m+1)*sizeof(int),&rows);
 53:   for (i=0; i<m; i++) {rows[i] = start + i;}

 55:   if (size == 1) {
 56:     MatCreateSeqDense(comm,M,M,PETSC_NULL,mat);
 57:   } else {
 58:     /*    MatCreateMPIDense(comm,m,M,M,M,PETSC_NULL,mat); */
 59:     MatCreateMPIAIJ(comm,m,m,M,M,0,0,0,0,mat);
 60:   }
 61: 
 62:   PCGetOperators(ksp->B,&A,PETSC_NULL,PETSC_NULL);

 64:   for (i=0; i<M; i++) {

 66:     VecSet(&zero,in);
 67:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
 68:     VecAssemblyBegin(in);
 69:     VecAssemblyEnd(in);

 71:     KSP_MatMult(ksp,A,in,out);
 72:     KSP_PCApply(ksp,ksp->B,out,in);
 73: 
 74:     VecGetArray(in,&array);
 75:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
 76:     VecRestoreArray(in,&array);

 78:   }
 79:   PetscFree(rows);
 80:   VecDestroy(in);
 81:   VecDestroy(out);
 82:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
 84:   return(0);
 85: }

 87:  #include petscblaslapack.h

 89: #undef __FUNCT__  
 91: /*@
 92:    KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the 
 93:    preconditioned operator using LAPACK.  

 95:    Collective on KSP

 97:    Input Parameter:
 98: +  ksp - iterative context obtained from KSPCreate()
 99: -  n - size of arrays r and c

101:    Output Parameters:
102: +  r - real part of computed eigenvalues
103: -  c - complex part of computed eigenvalues

105:    Notes:
106:    This approach is very slow but will generally provide accurate eigenvalue
107:    estimates.  This routine explicitly forms a dense matrix representing 
108:    the preconditioned operator, and thus will run only for relatively small
109:    problems, say n < 500.

111:    Many users may just want to use the monitoring routine
112:    KSPSingularValueMonitor() (which can be set with option -ksp_singmonitor)
113:    to print the singular values at each iteration of the linear solve.

115:    Level: advanced

117: .keywords: KSP, compute, eigenvalues, explicitly

119: .seealso: KSPComputeEigenvalues(), KSPSingularValueMonitor(), KSPComputeExtremeSingularValues()
120: @*/
121: int KSPComputeEigenvaluesExplicitly(KSP ksp,int nmax,PetscReal *r,PetscReal *c)
122: {
123:   Mat          BA;
124:   int          i,n,ierr,size,rank,dummy;
125:   MPI_Comm     comm = ksp->comm;
126:   PetscScalar  *array;
127:   Mat          A;
128:   int          m,row,nz,*cols;
129:   PetscScalar  *vals;

132:    KSPComputeExplicitOperator(ksp,&BA);
133:   MPI_Comm_size(comm,&size);
134:   MPI_Comm_rank(comm,&rank);

136:   ierr     = MatGetSize(BA,&n,&n);
137:   if (size > 1) { /* assemble matrix on first processor */
138:     if (!rank) {
139:       MatCreateMPIDense(ksp->comm,n,n,n,n,PETSC_NULL,&A);
140:     }
141:     else {
142:       MatCreateMPIDense(ksp->comm,0,n,n,n,PETSC_NULL,&A);
143:     }
144:     PetscLogObjectParent(BA,A);

146:     MatGetOwnershipRange(BA,&row,&dummy);
147:     MatGetLocalSize(BA,&m,&dummy);
148:     for (i=0; i<m; i++) {
149:       MatGetRow(BA,row,&nz,&cols,&vals);
150:       MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
151:       MatRestoreRow(BA,row,&nz,&cols,&vals);
152:       row++;
153:     }

155:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
156:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
157:     MatGetArray(A,&array);
158:   } else {
159:     MatGetArray(BA,&array);
160:   }

162: #if defined(PETSC_HAVE_ESSL)
163:   /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
164:   if (!rank) {
165:     PetscScalar    sdummy,*cwork;
166:     PetscReal *work,*realpart;
167:     int       clen,idummy,lwork,*perm,zero;

169: #if !defined(PETSC_USE_COMPLEX)
170:     clen = n;
171: #else
172:     clen = 2*n;
173: #endif
174:     ierr   = PetscMalloc(clen*sizeof(PetscScalar),&cwork);
175:     idummy = n;
176:     lwork  = 5*n;
177:     ierr   = PetscMalloc(lwork*sizeof(PetscReal),&work);
178:     ierr   = PetscMalloc(n*sizeof(PetscReal),&realpart);
179:     zero   = 0;
180:     LAgeev_(&zero,array,&n,cwork,&sdummy,&idummy,&idummy,&n,work,&lwork);
181:     PetscFree(work);

183:     /* For now we stick with the convention of storing the real and imaginary
184:        components of evalues separately.  But is this what we really want? */
185:     PetscMalloc(n*sizeof(int),&perm);

187: #if !defined(PETSC_USE_COMPLEX)
188:     for (i=0; i<n; i++) {
189:       realpart[i] = cwork[2*i];
190:       perm[i]     = i;
191:     }
192:     PetscSortRealWithPermutation(n,realpart,perm);
193:     for (i=0; i<n; i++) {
194:       r[i] = cwork[2*perm[i]];
195:       c[i] = cwork[2*perm[i]+1];
196:     }
197: #else
198:     for (i=0; i<n; i++) {
199:       realpart[i] = PetscRealPart(cwork[i]);
200:       perm[i]     = i;
201:     }
202:     PetscSortRealWithPermutation(n,realpart,perm);
203:     for (i=0; i<n; i++) {
204:       r[i] = PetscRealPart(cwork[perm[i]]);
205:       c[i] = PetscImaginaryPart(cwork[perm[i]]);
206:     }
207: #endif
208:     PetscFree(perm);
209:     PetscFree(realpart);
210:     PetscFree(cwork);
211:   }
212: #elif !defined(PETSC_USE_COMPLEX)
213:   if (!rank) {
214:     PetscScalar    *work,sdummy;
215:     PetscReal *realpart,*imagpart;
216:     int       idummy,lwork,*perm;

218:     idummy   = n;
219:     lwork    = 5*n;
220:     ierr     = PetscMalloc(2*n*sizeof(PetscReal),&realpart);
221:     imagpart = realpart + n;
222:     ierr     = PetscMalloc(5*n*sizeof(PetscReal),&work);
223: #if defined(PETSC_MISSING_LAPACK_GEEV) 
224:   SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavilablenNot able to provide eigen values.");
225: #else
226:     LAgeev_("N","N",&n,array,&n,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&ierr);
227: #endif
228:     if (ierr) SETERRQ1(PETSC_ERR_LIB,"Error in LAPACK routine %d",ierr);
229:     PetscFree(work);
230:     PetscMalloc(n*sizeof(int),&perm);
231:     for (i=0; i<n; i++) { perm[i] = i;}
232:     PetscSortRealWithPermutation(n,realpart,perm);
233:     for (i=0; i<n; i++) {
234:       r[i] = realpart[perm[i]];
235:       c[i] = imagpart[perm[i]];
236:     }
237:     PetscFree(perm);
238:     PetscFree(realpart);
239:   }
240: #else
241:   if (!rank) {
242:     PetscScalar *work,sdummy,*eigs;
243:     PetscReal *rwork;
244:     int    idummy,lwork,*perm;

246:     idummy   = n;
247:     lwork    = 5*n;
248:     PetscMalloc(5*n*sizeof(PetscScalar),&work);
249:     PetscMalloc(2*n*sizeof(PetscReal),&rwork);
250:     PetscMalloc(n*sizeof(PetscScalar),&eigs);
251: #if defined(PETSC_MISSING_LAPACK_GEEV) 
252:   SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavilablenNot able to provide eigen values.");
253: #else
254:     LAgeev_("N","N",&n,array,&n,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&ierr);
255: #endif
256:     if (ierr) SETERRQ1(PETSC_ERR_LIB,"Error in LAPACK routine %d",ierr);
257:     PetscFree(work);
258:     PetscFree(rwork);
259:     PetscMalloc(n*sizeof(int),&perm);
260:     for (i=0; i<n; i++) { perm[i] = i;}
261:     for (i=0; i<n; i++) { r[i]    = PetscRealPart(eigs[i]);}
262:     PetscSortRealWithPermutation(n,r,perm);
263:     for (i=0; i<n; i++) {
264:       r[i] = PetscRealPart(eigs[perm[i]]);
265:       c[i] = PetscImaginaryPart(eigs[perm[i]]);
266:     }
267:     PetscFree(perm);
268:     PetscFree(eigs);
269:   }
270: #endif  
271:   if (size > 1) {
272:     MatRestoreArray(A,&array);
273:     MatDestroy(A);
274:   } else {
275:     MatRestoreArray(BA,&array);
276:   }
277:   MatDestroy(BA);
278:   return(0);
279: }