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 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 26 ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr); 27 ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr); 28 ierr = PetscBLASIntCast(N,&idummy);CHKERRQ(ierr); 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 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 75 ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr); 76 idummy = -1; /* unused */ 77 ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr); 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 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 137 ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr); 138 ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr); 139 ierr = PetscBLASIntCast(N,&idummy);CHKERRQ(ierr); 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) PetscFunctionReturn(0); 144 145 /* copy R matrix to work space */ 146 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 147 148 /* compute eigenvalues */ 149 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 150 LAPACKgeev_("N","N",&bn,R,&bN,realpart,imagpart,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,&lierr); 151 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine %d",(int)lierr); 152 ierr = PetscFPTrapPop();CHKERRQ(ierr); 153 ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr); 154 for (i=0; i<n; i++) { perm[i] = i;} 155 ierr = PetscSortRealWithPermutation(n,realpart,perm);CHKERRQ(ierr); 156 for (i=0; i<n; i++) { 157 r[i] = realpart[perm[i]]; 158 c[i] = imagpart[perm[i]]; 159 } 160 ierr = PetscFree(perm);CHKERRQ(ierr); 161 #else 162 KSP_GMRES *gmres = (KSP_GMRES*)ksp->data; 163 PetscErrorCode ierr; 164 PetscInt n = gmres->it + 1,N = gmres->max_k + 1,i,*perm; 165 PetscScalar *R = gmres->Rsvd,*work = R + N*N,*eigs = work + 5*N,sdummy; 166 PetscBLASInt bn,bN,lwork,idummy,lierr; 167 168 PetscFunctionBegin; 169 ierr = PetscBLASIntCast(n,&bn);CHKERRQ(ierr); 170 ierr = PetscBLASIntCast(N,&bN);CHKERRQ(ierr); 171 ierr = PetscBLASIntCast(5*N,&lwork);CHKERRQ(ierr); 172 ierr = PetscBLASIntCast(N,&idummy);CHKERRQ(ierr); 173 if (nmax < n) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_SIZ,"Not enough room in work space r and c for eigenvalues"); 174 *neig = n; 175 176 if (!n) { 177 PetscFunctionReturn(0); 178 } 179 /* copy R matrix to work space */ 180 ierr = PetscMemcpy(R,gmres->hes_origin,N*N*sizeof(PetscScalar));CHKERRQ(ierr); 181 182 /* compute eigenvalues */ 183 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 184 LAPACKgeev_("N","N",&bn,R,&bN,eigs,&sdummy,&idummy,&sdummy,&idummy,work,&lwork,gmres->Dsvd,&lierr); 185 if (lierr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in LAPACK routine"); 186 ierr = PetscFPTrapPop();CHKERRQ(ierr); 187 ierr = PetscMalloc(n*sizeof(PetscInt),&perm);CHKERRQ(ierr); 188 for (i=0; i<n; i++) { perm[i] = i;} 189 for (i=0; i<n; i++) { r[i] = PetscRealPart(eigs[i]);} 190 ierr = PetscSortRealWithPermutation(n,r,perm);CHKERRQ(ierr); 191 for (i=0; i<n; i++) { 192 r[i] = PetscRealPart(eigs[perm[i]]); 193 c[i] = PetscImaginaryPart(eigs[perm[i]]); 194 } 195 ierr = PetscFree(perm);CHKERRQ(ierr); 196 #endif 197 PetscFunctionReturn(0); 198 } 199 200 201 202 203