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