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, bN ,lwork, idummy,lierr; 22 PetscScalar *R = gmres->Rsvd,*work = R + N*N,sdummy; 23 PetscReal *realpart = gmres->Dsvd; 24 25 PetscFunctionBegin; 26 bn = PetscBLASIntCast(n); 27 bN = PetscBLASIntCast(N); 28 lwork = PetscBLASIntCast(5*N); 29 idummy = PetscBLASIntCast(N); 30 if (!n) { 31 *emax = *emin = 1.0; 32 PetscFunctionReturn(0); 33 } 34 /* copy R matrix to work space */ 35 ierr = PetscMemcpy(R,gmres->hh_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 36 37 /* zero below diagonal garbage */ 38 for (i=0; i<n; i++) { 39 R[i*N+i+1] = 0.0; 40 } 41 42 /* compute Singular Values */ 43 #if !defined(PETSC_USE_COMPLEX) 44 LAPACKgesvd_("N","N",&bn,&bn,R,&bN,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr); 45 #else 46 LAPACKgesvd_("N","N",&bn,&bn,R,&bN,realpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,realpart+N,&lierr); 47 #endif 48 if (lierr) SETERRQ1(PETSC_ERR_LIB,"Error in SVD Lapack routine %d",(int)lierr); 49 50 *emin = realpart[n-1]; 51 *emax = realpart[0]; 52 #endif 53 PetscFunctionReturn(0); 54 } 55 56 /* ------------------------------------------------------------------------ */ 57 /* ESSL has a different calling sequence for dgeev() and zgeev() than standard LAPACK */ 58 #undef __FUNCT__ 59 #define __FUNCT__ "KSPComputeEigenvalues_GMRES" 60 PetscErrorCode KSPComputeEigenvalues_GMRES(KSP ksp,PetscInt nmax,PetscReal *r,PetscReal *c,PetscInt *neig) 61 { 62 #if defined(PETSC_HAVE_ESSL) 63 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 64 PetscErrorCode ierr; 65 PetscInt n = gmres->it + 1,N = gmres->max_k + 1,lwork = 5*N; 66 PetscInt i,*perm; 67 PetscScalar *R = gmres->Rsvd; 68 PetscScalar *cwork = R + N*N,sdummy; 69 PetscReal *work,*realpart = gmres->Dsvd ; 70 PetscBLASInt zero = 0,idummy = PetscBLASIntCast(N); 71 72 PetscFunctionBegin; 73 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 74 *neig = n; 75 76 if (!n) { 77 PetscFunctionReturn(0); 78 } 79 /* copy R matrix to work space */ 80 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 81 82 /* compute eigenvalues */ 83 84 /* for ESSL version need really cwork of length N (complex), 2N 85 (real); already at least 5N of space has been allocated */ 86 87 ierr = PetscMalloc(lwork*sizeof(PetscReal),&work);CHKERRQ(ierr); 88 LAPACKgeev_(&zero,R,&idummy,cwork,&sdummy,&idummy,&idummy,&n,work,&lwork); 89 ierr = PetscFree(work);CHKERRQ(ierr); 90 91 /* For now we stick with the convention of storing the real and imaginary 92 components of evalues separately. But is this what we really want? */ 93 ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr); 94 95 #if !defined(PETSC_USE_COMPLEX) 96 for (i=0; i<n; i++) { 97 realpart[i] = cwork[2*i]; 98 perm[i] = i; 99 } 100 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 101 for (i=0; i<n; i++) { 102 r[i] = cwork[2*perm[i]]; 103 c[i] = cwork[2*perm[i]+1]; 104 } 105 #else 106 for (i=0; i<n; i++) { 107 realpart[i] = PetscRealPart(cwork[i]); 108 perm[i] = i; 109 } 110 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 111 for (i=0; i<n; i++) { 112 r[i] = PetscRealPart(cwork[perm[i]]); 113 c[i] = PetscImaginaryPart(cwork[perm[i]]); 114 } 115 #endif 116 ierr = PetscFree(perm);CHKERRQ(ierr); 117 #elif defined(PETSC_MISSING_LAPACK_GEEV) 118 PetscFunctionBegin; 119 SETERRQ(PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values."); 120 #elif !defined(PETSC_USE_COMPLEX) 121 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 122 PetscErrorCode ierr; 123 PetscInt n = gmres->it + 1,N = gmres->max_k + 1,i,*perm; 124 PetscBLASInt bn, bN, lwork, idummy, lierr; 125 PetscScalar *R = gmres->Rsvd,*work = R + N*N; 126 PetscScalar *realpart = gmres->Dsvd,*imagpart = realpart + N,sdummy; 127 128 PetscFunctionBegin; 129 bn = PetscBLASIntCast(n); 130 bN = PetscBLASIntCast(N); 131 lwork = PetscBLASIntCast(5*N); 132 idummy = PetscBLASIntCast(N); 133 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 134 *neig = n; 135 136 if (!n) { 137 PetscFunctionReturn(0); 138 } 139 140 /* copy R matrix to work space */ 141 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 142 143 /* compute eigenvalues */ 144 LAPACKgeev_("N","N",&bn,R,&bN,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr); 145 if (lierr) SETERRQ1(PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr); 146 ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr); 147 for (i=0; i<n; i++) { perm[i] = i;} 148 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 149 for (i=0; i<n; i++) { 150 r[i] = realpart[perm[i]]; 151 c[i] = imagpart[perm[i]]; 152 } 153 ierr = PetscFree(perm);CHKERRQ(ierr); 154 #else 155 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 156 PetscErrorCode ierr; 157 PetscInt n = gmres->it + 1,N = gmres->max_k + 1,i,*perm; 158 PetscScalar *R = gmres->Rsvd,*work = R + N*N,*eigs = work + 5*N,sdummy; 159 PetscBLASInt bn,bN,lwork,idummy,lierr; 160 161 PetscFunctionBegin; 162 bn = PetscBLASIntCast(n); 163 bN = PetscBLASIntCast(N); 164 lwork = PetscBLASIntCast(5*N); 165 idummy = PetscBLASIntCast(N); 166 if (nmax < n) SETERRQ(PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 167 *neig = n; 168 169 if (!n) { 170 PetscFunctionReturn(0); 171 } 172 /* copy R matrix to work space */ 173 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 174 175 /* compute eigenvalues */ 176 LAPACKgeev_("N","N",&bn,R,&bN,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,gmres->Dsvd,&lierr); 177 if (lierr) SETERRQ(PETSC_ERR_LIB,"Error in LAPACK routine"); 178 ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr); 179 for (i=0; i<n; i++) { perm[i] = i;} 180 for (i=0; i<n; i++) { r[i] = PetscRealPart(eigs[i]);} 181 ierr = PetscSortRealWithPermutation(n,r,perm);CHKERRQ(ierr); 182 for (i=0; i<n; i++) { 183 r[i] = PetscRealPart(eigs[perm[i]]); 184 c[i] = PetscImaginaryPart(eigs[perm[i]]); 185 } 186 ierr = PetscFree(perm);CHKERRQ(ierr); 187 #endif 188 PetscFunctionReturn(0); 189 } 190 191 192 193 194