Actual source code: eige.c

      1: #define PETSCKSP_DLL

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

      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: PetscErrorCode PETSCKSP_DLLEXPORT KSPComputeExplicitOperator(KSP ksp,Mat *mat)
     33: {
     34:   Vec            in,out;
     36:   PetscMPIInt    size;
     37:   PetscInt       i,M,m,*rows,start,end;
     38:   Mat            A;
     39:   MPI_Comm       comm;
     40:   PetscScalar    *array,zero = 0.0,one = 1.0;

     45:   comm = ksp->comm;

     47:   MPI_Comm_size(comm,&size);

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

     57:   MatCreate(comm,mat);
     58:   MatSetSizes(*mat,m,m,M,M);
     59:   if (size == 1) {
     60:     MatSetType(*mat,MATSEQDENSE);
     61:     MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
     62:   } else {
     63:     MatSetType(*mat,MATMPIAIJ);
     64:     MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);
     65:   }
     66:
     67:   PCGetOperators(ksp->pc,&A,PETSC_NULL,PETSC_NULL);

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

     71:     VecSet(in,zero);
     72:     VecSetValues(in,1,&i,&one,INSERT_VALUES);
     73:     VecAssemblyBegin(in);
     74:     VecAssemblyEnd(in);

     76:     KSP_MatMult(ksp,A,in,out);
     77:     KSP_PCApply(ksp,out,in);
     78:
     79:     VecGetArray(in,&array);
     80:     MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
     81:     VecRestoreArray(in,&array);

     83:   }
     84:   PetscFree(rows);
     85:   VecDestroy(in);
     86:   VecDestroy(out);
     87:   MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
     88:   MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
     89:   return(0);
     90: }

     92:  #include petscblaslapack.h

     96: /*@
     97:    KSPComputeEigenvaluesExplicitly - Computes all of the eigenvalues of the
     98:    preconditioned operator using LAPACK.

    100:    Collective on KSP

    102:    Input Parameter:
    103: +  ksp - iterative context obtained from KSPCreate()
    104: -  n - size of arrays r and c

    106:    Output Parameters:
    107: +  r - real part of computed eigenvalues
    108: -  c - complex part of computed eigenvalues

    110:    Notes:
    111:    This approach is very slow but will generally provide accurate eigenvalue
    112:    estimates.  This routine explicitly forms a dense matrix representing
    113:    the preconditioned operator, and thus will run only for relatively small
    114:    problems, say n < 500.

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

    120:    The preconditoner operator, rhs vector, solution vectors should be
    121:    set before this routine is called. i.e use KSPSetOperators(),KSPSolve() or
    122:    KSPSetOperators()

    124:    Level: advanced

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

    128: .seealso: KSPComputeEigenvalues(), KSPSingularValueMonitor(), KSPComputeExtremeSingularValues(), KSPSetOperators(), KSPSolve()
    129: @*/
    130: PetscErrorCode PETSCKSP_DLLEXPORT KSPComputeEigenvaluesExplicitly(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c)
    131: {
    132:   Mat                BA;
    133:   PetscErrorCode     ierr;
    134:   PetscMPIInt        size,rank;
    135:   MPI_Comm           comm = ksp->comm;
    136:   PetscScalar        *array;
    137:   Mat                A;
    138:   PetscInt           m,row,nz,i,n,dummy;
    139:   const PetscInt     *cols;
    140:   const PetscScalar  *vals;

    143:   KSPComputeExplicitOperator(ksp,&BA);
    144:   MPI_Comm_size(comm,&size);
    145:   MPI_Comm_rank(comm,&rank);

    147:   MatGetSize(BA,&n,&n);
    148:   if (size > 1) { /* assemble matrix on first processor */
    149:     MatCreate(ksp->comm,&A);
    150:     if (!rank) {
    151:       MatSetSizes(A,n,n,n,n);
    152:     } else {
    153:       MatSetSizes(A,0,0,n,n);
    154:     }
    155:     MatSetType(A,MATMPIDENSE);
    156:     MatMPIDenseSetPreallocation(A,PETSC_NULL);
    157:     PetscLogObjectParent(BA,A);

    159:     MatGetOwnershipRange(BA,&row,&dummy);
    160:     MatGetLocalSize(BA,&m,&dummy);
    161:     for (i=0; i<m; i++) {
    162:       MatGetRow(BA,row,&nz,&cols,&vals);
    163:       MatSetValues(A,1,&row,nz,cols,vals,INSERT_VALUES);
    164:       MatRestoreRow(BA,row,&nz,&cols,&vals);
    165:       row++;
    166:     }

    168:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    169:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    170:     MatGetArray(A,&array);
    171:   } else {
    172:     MatGetArray(BA,&array);
    173:   }

    175: #if defined(PETSC_HAVE_ESSL)
    176:   /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
    177:   if (!rank) {
    178:     PetscScalar sdummy,*cwork;
    179:     PetscReal   *work,*realpart;
    180:     PetscInt         clen,idummy,lwork,*perm,zero;

    182: #if !defined(PETSC_USE_COMPLEX)
    183:     clen = n;
    184: #else
    185:     clen = 2*n;
    186: #endif
    187:     PetscMalloc(clen*sizeof(PetscScalar),&cwork);
    188:     idummy = n;
    189:     lwork  = 5*n;
    190:     PetscMalloc(lwork*sizeof(PetscReal),&work);
    191:     PetscMalloc(n*sizeof(PetscReal),&realpart);
    192:     zero   = 0;
    193:     LAPACKgeev_(&zero,array,&n,cwork,&sdummy,&idummy,&idummy,&n,work,&lwork);
    194:     PetscFree(work);

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

    200: #if !defined(PETSC_USE_COMPLEX)
    201:     for (i=0; i<n; i++) {
    202:       realpart[i] = cwork[2*i];
    203:       perm[i]     = i;
    204:     }
    205:     PetscSortRealWithPermutation(n,realpart,perm);
    206:     for (i=0; i<n; i++) {
    207:       r[i] = cwork[2*perm[i]];
    208:       c[i] = cwork[2*perm[i]+1];
    209:     }
    210: #else
    211:     for (i=0; i<n; i++) {
    212:       realpart[i] = PetscRealPart(cwork[i]);
    213:       perm[i]     = i;
    214:     }
    215:     PetscSortRealWithPermutation(n,realpart,perm);
    216:     for (i=0; i<n; i++) {
    217:       r[i] = PetscRealPart(cwork[perm[i]]);
    218:       c[i] = PetscImaginaryPart(cwork[perm[i]]);
    219:     }
    220: #endif
    221:     PetscFree(perm);
    222:     PetscFree(realpart);
    223:     PetscFree(cwork);
    224:   }
    225: #elif !defined(PETSC_USE_COMPLEX)
    226:   if (!rank) {
    227:     PetscScalar  *work;
    228:     PetscReal    *realpart,*imagpart;
    229:     PetscBLASInt idummy,lwork;
    230:     PetscInt     *perm;

    232:     idummy   = n;
    233:     lwork    = 5*n;
    234:     PetscMalloc(2*n*sizeof(PetscReal),&realpart);
    235:     imagpart = realpart + n;
    236:     PetscMalloc(5*n*sizeof(PetscReal),&work);
    237: #if defined(PETSC_MISSING_LAPACK_GEEV)
    238:     SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
    239: #else
    240:     {
    241:       PetscBLASInt lierr;
    242:       PetscBLASInt bn = (PetscBLASInt) n;
    243:       PetscScalar sdummy;
    244:       LAPACKgeev_("N","N",&bn,array,&bn,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr);
    245:       if (lierr) SETERRQ1(PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
    246:     }
    247: #endif
    248:     PetscFree(work);
    249:     PetscMalloc(n*sizeof(PetscInt),&perm);
    250:     for (i=0; i<n; i++) { perm[i] = i;}
    251:     PetscSortRealWithPermutation(n,realpart,perm);
    252:     for (i=0; i<n; i++) {
    253:       r[i] = realpart[perm[i]];
    254:       c[i] = imagpart[perm[i]];
    255:     }
    256:     PetscFree(perm);
    257:     PetscFree(realpart);
    258:   }
    259: #else
    260:   if (!rank) {
    261:     PetscScalar  *work,*eigs;
    262:     PetscReal    *rwork;
    263:     PetscBLASInt idummy,lwork;
    264:     PetscInt     *perm;

    266:     idummy   = n;
    267:     lwork    = 5*n;
    268:     PetscMalloc(5*n*sizeof(PetscScalar),&work);
    269:     PetscMalloc(2*n*sizeof(PetscReal),&rwork);
    270:     PetscMalloc(n*sizeof(PetscScalar),&eigs);
    271: #if defined(PETSC_MISSING_LAPACK_GEEV)
    272:     SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
    273: #else
    274:     {
    275:       PetscBLASInt lierr;
    276:       PetscScalar sdummy;
    277:       LAPACKgeev_("N","N",&n,array,&n,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,rwork,&lierr);
    278:       if (lierr) SETERRQ1(PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
    279:     }
    280: #endif
    281:     PetscFree(work);
    282:     PetscFree(rwork);
    283:     PetscMalloc(n*sizeof(PetscInt),&perm);
    284:     for (i=0; i<n; i++) { perm[i] = i;}
    285:     for (i=0; i<n; i++) { r[i]    = PetscRealPart(eigs[i]);}
    286:     PetscSortRealWithPermutation(n,r,perm);
    287:     for (i=0; i<n; i++) {
    288:       r[i] = PetscRealPart(eigs[perm[i]]);
    289:       c[i] = PetscImaginaryPart(eigs[perm[i]]);
    290:     }
    291:     PetscFree(perm);
    292:     PetscFree(eigs);
    293:   }
    294: #endif
    295:   if (size > 1) {
    296:     MatRestoreArray(A,&array);
    297:     MatDestroy(A);
    298:   } else {
    299:     MatRestoreArray(BA,&array);
    300:   }
    301:   MatDestroy(BA);
    302:   return(0);
    303: }