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