1 #define PETSCKSP_DLL 2 3 #include "src/ksp/ksp/impls/gmres/gmresp.h" 4 #include "petscblaslapack.h" 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "KSPComputeExtremeSingularValues_GMRES" 8 PetscErrorCode KSPComputeExtremeSingularValues_GMRES(KSP ksp,PetscReal *emax,PetscReal *emin) 9 { 10 #if defined(PETSC_MISSING_LAPACK_GESVD) 11 PetscFunctionBegin; 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; 19 PetscErrorCode ierr; 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; 24 25 PetscFunctionBegin; 26 if (!n) { 27 *emax = *emin = 1.0; 28 PetscFunctionReturn(0); 29 } 30 /* copy R matrix to work space */ 31 ierr = PetscMemcpy(R,gmres->hh_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 32 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); 45 46 *emin = realpart[n-1]; 47 *emax = realpart[0]; 48 #endif 49 PetscFunctionReturn(0); 50 } 51 52 /* ------------------------------------------------------------------------ */ 53 /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */ 54 #undef __FUNCT__ 55 #define __FUNCT__ "KSPComputeEigenvalues_GMRES" 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; 60 PetscErrorCode ierr; 61 PetscInt n = gmres->it + 1,N = gmres->max_k + 1,lwork = 5*N; 62 PetscInt i,*perm; 63 PetscBLASInt zero = 0,idummy = N; 64 PetscScalar *R = gmres->Rsvd; 65 PetscScalar *cwork = R + N*N,sdummy; 66 PetscReal *work,*realpart = gmres->Dsvd ; 67 68 PetscFunctionBegin; 69 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 70 *neig = n; 71 72 if (!n) { 73 PetscFunctionReturn(0); 74 } 75 /* copy R matrix to work space */ 76 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 77 78 /* compute eigenvalues */ 79 80 /* for ESSL version need really cwork of length N (complex), 2N 81 (real); already at least 5N of space has been allocated */ 82 83 ierr = PetscMalloc(lwork*sizeof(PetscReal),&work);CHKERRQ(ierr); 84 LAPACKgeev_(&zero,R,&idummy,cwork,&sdummy,&idummy,&idummy,&n,work,&lwork); 85 ierr = PetscFree(work);CHKERRQ(ierr); 86 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 ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr); 90 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 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 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 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 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 ierr = PetscFree(perm);CHKERRQ(ierr); 113 #elif defined(PETSC_MISSING_LAPACK_GEEV) 114 PetscFunctionBegin; 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; 118 PetscErrorCode ierr; 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; 123 124 PetscFunctionBegin; 125 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 126 *neig = n; 127 128 if (!n) { 129 PetscFunctionReturn(0); 130 } 131 132 /* copy R matrix to work space */ 133 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 134 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 ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr); 139 for (i=0; i<n; i++) { perm[i] = i;} 140 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 141 for (i=0; i<n; i++) { 142 r[i] = realpart[perm[i]]; 143 c[i] = imagpart[perm[i]]; 144 } 145 ierr = PetscFree(perm);CHKERRQ(ierr); 146 #else 147 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 148 PetscErrorCode ierr; 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; 151 152 PetscFunctionBegin; 153 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 154 *neig = n; 155 156 if (!n) { 157 PetscFunctionReturn(0); 158 } 159 /* copy R matrix to work space */ 160 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 161 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 ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr); 166 for (i=0; i<n; i++) { perm[i] = i;} 167 for (i=0; i<n; i++) { r[i] = PetscRealPart(eigs[i]);} 168 ierr = PetscSortRealWithPermutation(n,r,perm);CHKERRQ(ierr); 169 for (i=0; i<n; i++) { 170 r[i] = PetscRealPart(eigs[perm[i]]); 171 c[i] = PetscImaginaryPart(eigs[perm[i]]); 172 } 173 ierr = PetscFree(perm);CHKERRQ(ierr); 174 #endif 175 PetscFunctionReturn(0); 176 } 177 178 179 180 181