1 #define PETSCKSP_DLL 2 3 #include "private/pcimpl.h" /*I "petscpc.h" I*/ 4 #include "petscblaslapack.h" 5 6 /* 7 Private context (data structure) for the SVD preconditioner. 8 */ 9 typedef struct { 10 Vec diag,work; 11 Mat A,U,V; 12 PetscInt nzero; 13 } PC_SVD; 14 15 16 /* -------------------------------------------------------------------------- */ 17 /* 18 PCSetUp_SVD - Prepares for the use of the SVD preconditioner 19 by setting data structures and options. 20 21 Input Parameter: 22 . pc - the preconditioner context 23 24 Application Interface Routine: PCSetUp() 25 26 Notes: 27 The interface routine PCSetUp() is not usually called directly by 28 the user, but instead is called by PCApply() if necessary. 29 */ 30 #undef __FUNCT__ 31 #define __FUNCT__ "PCSetUp_SVD" 32 static PetscErrorCode PCSetUp_SVD(PC pc) 33 { 34 #if defined(PETSC_MISSING_LAPACK_GESVD) 35 SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable\nNot able to provide singular value estimates."); 36 #else 37 PC_SVD *jac = (PC_SVD*)pc->data; 38 PetscErrorCode ierr; 39 PetscScalar *a,*u,*v,*d,*work; 40 PetscBLASInt nb,lwork,lierr; 41 PetscInt i,n; 42 43 PetscFunctionBegin; 44 if (!jac->diag) { 45 /* assume square matrices */ 46 ierr = MatGetVecs(pc->mat,&jac->diag,&jac->work);CHKERRQ(ierr); 47 } 48 if (jac->A) { 49 ierr = MatDestroy(jac->A);CHKERRQ(ierr); 50 } 51 ierr = MatConvert(pc->mat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);CHKERRQ(ierr); 52 if (!jac->U) { 53 ierr = MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);CHKERRQ(ierr); 54 ierr = MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->V);CHKERRQ(ierr); 55 } 56 ierr = MatGetSize(pc->mat,&n,PETSC_NULL);CHKERRQ(ierr); 57 nb = PetscBLASIntCast(n); 58 lwork = 5*nb; 59 ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr); 60 ierr = MatGetArray(jac->A,&a);CHKERRQ(ierr); 61 ierr = MatGetArray(jac->U,&u);CHKERRQ(ierr); 62 ierr = MatGetArray(jac->V,&v);CHKERRQ(ierr); 63 ierr = VecGetArray(jac->diag,&d);CHKERRQ(ierr); 64 #if !defined(PETSC_USE_COMPLEX) 65 LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr); 66 #else 67 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex"); 68 #endif 69 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr); 70 ierr = MatRestoreArray(jac->A,&a);CHKERRQ(ierr); 71 ierr = MatRestoreArray(jac->U,&u);CHKERRQ(ierr); 72 ierr = MatRestoreArray(jac->V,&v);CHKERRQ(ierr); 73 jac->nzero = 0; 74 for (i=0; i<n; i++) { 75 if (PetscRealPart(d[i]) < 1.e-12) {jac->nzero = n - i;break;} 76 d[i] = 1.0/d[i]; 77 } 78 ierr = PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero); 79 ierr = VecRestoreArray(jac->diag,&d);CHKERRQ(ierr); 80 #if defined(foo) 81 { 82 PetscViewer viewer; 83 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr); 84 ierr = MatView(jac->A,viewer);CHKERRQ(ierr); 85 ierr = MatView(jac->U,viewer);CHKERRQ(ierr); 86 ierr = MatView(jac->V,viewer);CHKERRQ(ierr); 87 ierr = VecView(jac->diag,viewer);CHKERRQ(ierr); 88 ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr); 89 } 90 #endif 91 ierr = PetscFree(work); 92 PetscFunctionReturn(0); 93 #endif 94 } 95 96 /* -------------------------------------------------------------------------- */ 97 /* 98 PCApply_SVD - Applies the SVD preconditioner to a vector. 99 100 Input Parameters: 101 . pc - the preconditioner context 102 . x - input vector 103 104 Output Parameter: 105 . y - output vector 106 107 Application Interface Routine: PCApply() 108 */ 109 #undef __FUNCT__ 110 #define __FUNCT__ "PCApply_SVD" 111 static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y) 112 { 113 PC_SVD *jac = (PC_SVD*)pc->data; 114 Vec work = jac->work; 115 PetscErrorCode ierr; 116 117 PetscFunctionBegin; 118 ierr = MatMultTranspose(jac->U,x,work);CHKERRQ(ierr); 119 ierr = VecPointwiseMult(work,work,jac->diag);CHKERRQ(ierr); 120 ierr = MatMultTranspose(jac->V,work,y);CHKERRQ(ierr); 121 PetscFunctionReturn(0); 122 } 123 124 /* -------------------------------------------------------------------------- */ 125 /* 126 PCDestroy_SVD - Destroys the private context for the SVD preconditioner 127 that was created with PCCreate_SVD(). 128 129 Input Parameter: 130 . pc - the preconditioner context 131 132 Application Interface Routine: PCDestroy() 133 */ 134 #undef __FUNCT__ 135 #define __FUNCT__ "PCDestroy_SVD" 136 static PetscErrorCode PCDestroy_SVD(PC pc) 137 { 138 PC_SVD *jac = (PC_SVD*)pc->data; 139 PetscErrorCode ierr; 140 141 PetscFunctionBegin; 142 if (jac->A) { 143 ierr = MatDestroy(jac->A);CHKERRQ(ierr); 144 } 145 if (jac->U) { 146 ierr = MatDestroy(jac->U);CHKERRQ(ierr); 147 } 148 if (jac->V) { 149 ierr = MatDestroy(jac->V);CHKERRQ(ierr); 150 } 151 if (jac->diag) { 152 ierr = VecDestroy(jac->diag);CHKERRQ(ierr); 153 } 154 ierr = PetscFree(jac);CHKERRQ(ierr); 155 PetscFunctionReturn(0); 156 } 157 158 #undef __FUNCT__ 159 #define __FUNCT__ "PCSetFromOptions_SVD" 160 static PetscErrorCode PCSetFromOptions_SVD(PC pc) 161 { 162 PetscErrorCode ierr; 163 164 PetscFunctionBegin; 165 ierr = PetscOptionsHead("SVD options");CHKERRQ(ierr); 166 ierr = PetscOptionsTail();CHKERRQ(ierr); 167 PetscFunctionReturn(0); 168 } 169 170 /* -------------------------------------------------------------------------- */ 171 /* 172 PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD, 173 and sets this as the private data within the generic preconditioning 174 context, PC, that was created within PCCreate(). 175 176 Input Parameter: 177 . pc - the preconditioner context 178 179 Application Interface Routine: PCCreate() 180 */ 181 182 /*MC 183 PCSVD - Use pseudo inverse defined by SVD of operator 184 185 Level: advanced 186 187 Concepts: SVD 188 189 Zero entries along the diagonal are replaced with the value 0.0 190 191 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 192 M*/ 193 194 EXTERN_C_BEGIN 195 #undef __FUNCT__ 196 #define __FUNCT__ "PCCreate_SVD" 197 PetscErrorCode PCCreate_SVD(PC pc) 198 { 199 PC_SVD *jac; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 /* 204 Creates the private data structure for this preconditioner and 205 attach it to the PC object. 206 */ 207 ierr = PetscNewLog(pc,PC_SVD,&jac);CHKERRQ(ierr); 208 pc->data = (void*)jac; 209 210 /* 211 Set the pointers for the functions that are provided above. 212 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 213 are called, they will automatically call these functions. Note we 214 choose not to provide a couple of these functions since they are 215 not needed. 216 */ 217 pc->ops->apply = PCApply_SVD; 218 pc->ops->applytranspose = PCApply_SVD; 219 pc->ops->setup = PCSetUp_SVD; 220 pc->ops->destroy = PCDestroy_SVD; 221 pc->ops->setfromoptions = PCSetFromOptions_SVD; 222 pc->ops->view = 0; 223 pc->ops->applyrichardson = 0; 224 PetscFunctionReturn(0); 225 } 226 EXTERN_C_END 227 228