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 PetscReal zerosing; /* measure of smallest singular value treated as nonzero */ 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->pmat,&jac->diag,&jac->work);CHKERRQ(ierr); 47 } 48 ierr = MatDestroy(&jac->A);CHKERRQ(ierr); 49 ierr = MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);CHKERRQ(ierr); 50 if (!jac->U) { 51 ierr = MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);CHKERRQ(ierr); 52 ierr = MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->V);CHKERRQ(ierr); 53 } 54 ierr = MatGetSize(pc->pmat,&n,PETSC_NULL);CHKERRQ(ierr); 55 nb = PetscBLASIntCast(n); 56 lwork = 5*nb; 57 ierr = PetscMalloc(lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr); 58 ierr = MatGetArray(jac->A,&a);CHKERRQ(ierr); 59 ierr = MatGetArray(jac->U,&u);CHKERRQ(ierr); 60 ierr = MatGetArray(jac->V,&v);CHKERRQ(ierr); 61 ierr = VecGetArray(jac->diag,&d);CHKERRQ(ierr); 62 #if !defined(PETSC_USE_COMPLEX) 63 LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr); 64 #else 65 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex"); 66 #endif 67 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr); 68 ierr = MatRestoreArray(jac->A,&a);CHKERRQ(ierr); 69 ierr = MatRestoreArray(jac->U,&u);CHKERRQ(ierr); 70 ierr = MatRestoreArray(jac->V,&v);CHKERRQ(ierr); 71 jac->nzero = 0; 72 ierr = PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1])); 73 for (i=0; i<n; i++) { 74 if (i == 0 && PetscRealPart(d[i]) == 0.0) {jac->nzero = n;break;} 75 if (PetscRealPart(d[i]) < jac->zerosing /* PetscRealPart(d[0]) */) {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 #undef __FUNCT__ 125 #define __FUNCT__ "PCReset_SVD" 126 static PetscErrorCode PCReset_SVD(PC pc) 127 { 128 PC_SVD *jac = (PC_SVD*)pc->data; 129 PetscErrorCode ierr; 130 131 PetscFunctionBegin; 132 ierr = MatDestroy(&jac->A);CHKERRQ(ierr); 133 ierr = MatDestroy(&jac->U);CHKERRQ(ierr); 134 ierr = MatDestroy(&jac->V);CHKERRQ(ierr); 135 ierr = VecDestroy(&jac->diag);CHKERRQ(ierr); 136 ierr = VecDestroy(&jac->work);CHKERRQ(ierr); 137 PetscFunctionReturn(0); 138 } 139 140 /* -------------------------------------------------------------------------- */ 141 /* 142 PCDestroy_SVD - Destroys the private context for the SVD preconditioner 143 that was created with PCCreate_SVD(). 144 145 Input Parameter: 146 . pc - the preconditioner context 147 148 Application Interface Routine: PCDestroy() 149 */ 150 #undef __FUNCT__ 151 #define __FUNCT__ "PCDestroy_SVD" 152 static PetscErrorCode PCDestroy_SVD(PC pc) 153 { 154 PetscErrorCode ierr; 155 156 PetscFunctionBegin; 157 ierr = PCReset_SVD(pc);CHKERRQ(ierr); 158 ierr = PetscFree(pc->data);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "PCSetFromOptions_SVD" 164 static PetscErrorCode PCSetFromOptions_SVD(PC pc) 165 { 166 PetscErrorCode ierr; 167 PC_SVD *jac = (PC_SVD*)pc->data; 168 169 PetscFunctionBegin; 170 ierr = PetscOptionsHead("SVD options");CHKERRQ(ierr); 171 ierr = PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,PETSC_NULL);CHKERRQ(ierr); 172 ierr = PetscOptionsTail();CHKERRQ(ierr); 173 PetscFunctionReturn(0); 174 } 175 176 /* -------------------------------------------------------------------------- */ 177 /* 178 PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD, 179 and sets this as the private data within the generic preconditioning 180 context, PC, that was created within PCCreate(). 181 182 Input Parameter: 183 . pc - the preconditioner context 184 185 Application Interface Routine: PCCreate() 186 */ 187 188 /*MC 189 PCSVD - Use pseudo inverse defined by SVD of operator 190 191 Level: advanced 192 193 Concepts: SVD 194 195 Options Database: 196 . -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero 197 198 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 199 M*/ 200 201 EXTERN_C_BEGIN 202 #undef __FUNCT__ 203 #define __FUNCT__ "PCCreate_SVD" 204 PetscErrorCode PCCreate_SVD(PC pc) 205 { 206 PC_SVD *jac; 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 /* 211 Creates the private data structure for this preconditioner and 212 attach it to the PC object. 213 */ 214 ierr = PetscNewLog(pc,PC_SVD,&jac);CHKERRQ(ierr); 215 jac->zerosing = 1.e-12; 216 pc->data = (void*)jac; 217 218 /* 219 Set the pointers for the functions that are provided above. 220 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.) 221 are called, they will automatically call these functions. Note we 222 choose not to provide a couple of these functions since they are 223 not needed. 224 */ 225 pc->ops->apply = PCApply_SVD; 226 pc->ops->applytranspose = PCApply_SVD; 227 pc->ops->setup = PCSetUp_SVD; 228 pc->ops->reset = PCReset_SVD; 229 pc->ops->destroy = PCDestroy_SVD; 230 pc->ops->setfromoptions = PCSetFromOptions_SVD; 231 pc->ops->view = 0; 232 pc->ops->applyrichardson = 0; 233 PetscFunctionReturn(0); 234 } 235 EXTERN_C_END 236 237