Actual source code: gmreig.c

      1: #define PETSCKSP_DLL

     3:  #include src/ksp/ksp/impls/gmres/gmresp.h
     4:  #include petscblaslapack.h

      8: PetscErrorCode KSPComputeExtremeSingularValues_GMRES(KSP ksp,PetscReal *emax,PetscReal *emin)
      9: {
     10: #if defined(PETSC_MISSING_LAPACK_GESVD)
     12:   /*
     13:       The Cray math libraries on T3D/T3E, and early versions of Intel Math Kernel Libraries (MKL)
     14:       for PCs do not seem to have the DGESVD() lapack routines
     15:   */
     16:   SETERRQ(PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable\nNot able to provide singular value estimates.");
     17: #else
     18:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
     20:   PetscInt       n = gmres->it + 1,i,N = gmres->max_k + 2;
     21:   PetscBLASInt   bn = (PetscBLASInt)n,bN = (PetscBLASInt)N,lwork = (PetscBLASInt)5*N,idummy = (PetscBLASInt)N,lierr;
     22:   PetscScalar    *R = gmres->Rsvd,*work = R + N*N,sdummy;
     23:   PetscReal      *realpart = gmres->Dsvd;

     26:   if (!n) {
     27:     *emax = *emin = 1.0;
     28:     return(0);
     29:   }
     30:   /* copy R matrix to work space */
     31:   PetscMemcpy(R,gmres->hh_origin,N*N*sizeof(PetscScalar));

     33:   /* zero below diagonal garbage */
     34:   for (i=0; i<n; i++) {
     35:     R[i*N+i+1] = 0.0;
     36:   }
     37:
     38:   /* compute Singular Values */
     39: #if !defined(PETSC_USE_COMPLEX)
     40:   LAPACKgesvd_("N","N",&bn,&bn,R,&bN,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr);
     41: #else
     42:   LAPACKgesvd_("N","N",&bn,&bn,R,&bN,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,realpart+N,&lierr);
     43: #endif
     44:   if (lierr) SETERRQ1(PETSC_ERR_LIB,"Error in SVD Lapack routine %d",(int)lierr);

     46:   *emin = realpart[n-1];
     47:   *emax = realpart[0];
     48: #endif
     49:   return(0);
     50: }

     52: /* ------------------------------------------------------------------------ */
     53: /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */
     56: PetscErrorCode KSPComputeEigenvalues_GMRES(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig)
     57: {
     58: #if defined(PETSC_HAVE_ESSL)
     59:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
     61:   PetscInt       n = gmres->it + 1,N = gmres->max_k + 1,lwork = 5*N;
     62:   PetscInt       idummy = N,i,*perm,zero;
     63:   PetscScalar    *R = gmres->Rsvd;
     64:   PetscScalar    *cwork = R + N*N,sdummy;
     65:   PetscReal      *work,*realpart = gmres->Dsvd ;

     68:   if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues");
     69:   *neig = n;

     71:   if (!n) {
     72:     return(0);
     73:   }
     74:   /* copy R matrix to work space */
     75:   PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));

     77:   /* compute eigenvalues */

     79:   /* for ESSL version need really cwork of length N (complex), 2N
     80:      (real); already at least 5N of space has been allocated */

     82:   PetscMalloc(lwork*sizeof(PetscReal),&work);
     83:   zero = 0;
     84:   LAPACKgeev_(&zero,R,&N,cwork,&sdummy,&idummy,&idummy,&n,work,&lwork);
     85:   PetscFree(work);

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

     91: #if !defined(PETSC_USE_COMPLEX)
     92:   for (i=0; i<n; i++) {
     93:     realpart[i] = cwork[2*i];
     94:     perm[i]     = i;
     95:   }
     96:   PetscSortRealWithPermutation(n,realpart,perm);
     97:   for (i=0; i<n; i++) {
     98:     r[i] = cwork[2*perm[i]];
     99:     c[i] = cwork[2*perm[i]+1];
    100:   }
    101: #else
    102:   for (i=0; i<n; i++) {
    103:     realpart[i] = PetscRealPart(cwork[i]);
    104:     perm[i]     = i;
    105:   }
    106:   PetscSortRealWithPermutation(n,realpart,perm);
    107:   for (i=0; i<n; i++) {
    108:     r[i] = PetscRealPart(cwork[perm[i]]);
    109:     c[i] = PetscImaginaryPart(cwork[perm[i]]);
    110:   }
    111: #endif
    112:   PetscFree(perm);
    113: #elif defined(PETSC_MISSING_LAPACK_GEEV)
    115:   SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
    116: #elif !defined(PETSC_USE_COMPLEX)
    117:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
    119:   PetscInt       n = gmres->it + 1,N = gmres->max_k + 1,i,*perm;
    120:   PetscBLASInt   bn = (PetscBLASInt)n,bN = (PetscBLASInt)N,lwork = (PetscBLASInt)5*N,idummy = (PetscBLASInt)N,lierr;
    121:   PetscScalar    *R = gmres->Rsvd,*work = R + N*N;
    122:   PetscScalar    *realpart = gmres->Dsvd,*imagpart = realpart + N,sdummy;

    125:   if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues");
    126:   *neig = n;

    128:   if (!n) {
    129:     return(0);
    130:   }

    132:   /* copy R matrix to work space */
    133:   PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));

    135:   /* compute eigenvalues */
    136:   LAPACKgeev_("N","N",&bn,R,&bN,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr);
    137:   if (lierr) SETERRQ1(PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr);
    138:   PetscMalloc(n*sizeof(PetscInt),&perm);
    139:   for (i=0; i<n; i++) { perm[i] = i;}
    140:   PetscSortRealWithPermutation(n,realpart,perm);
    141:   for (i=0; i<n; i++) {
    142:     r[i] = realpart[perm[i]];
    143:     c[i] = imagpart[perm[i]];
    144:   }
    145:   PetscFree(perm);
    146: #else
    147:   KSP_GMRES      *gmres = (KSP_GMRES*)ksp->data;
    149:   PetscInt       n = gmres->it + 1,N = gmres->max_k + 1,lwork = 5*N,idummy = N,i,*perm;
    150:   PetscScalar    *R = gmres->Rsvd,*work = R + N*N,*eigs = work + 5*N,sdummy;

    153:   if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues");
    154:   *neig = n;

    156:   if (!n) {
    157:     return(0);
    158:   }
    159:   /* copy R matrix to work space */
    160:   PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));

    162:   /* compute eigenvalues */
    163:   LAPACKgeev_("N","N",&n,R,&N,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,gmres->Dsvd,&ierr);
    164:   if (ierr) SETERRQ(PETSC_ERR_LIB,"Error in LAPACK routine");
    165:   PetscMalloc(n*sizeof(PetscInt),&perm);
    166:   for (i=0; i<n; i++) { perm[i] = i;}
    167:   for (i=0; i<n; i++) { r[i]    = PetscRealPart(eigs[i]);}
    168:   PetscSortRealWithPermutation(n,r,perm);
    169:   for (i=0; i<n; i++) {
    170:     r[i] = PetscRealPart(eigs[perm[i]]);
    171:     c[i] = PetscImaginaryPart(eigs[perm[i]]);
    172:   }
    173:   PetscFree(perm);
    174: #endif
    175:   return(0);
    176: }