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