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