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