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