1 2 /* 3 Routines to project vectors out of null spaces. 4 */ 5 6 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 7 8 PetscClassId MAT_NULLSPACE_CLASSID; 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "MatNullSpaceSetFunction" 12 /*@C 13 MatNullSpaceSetFunction - set a function that removes a null space from a vector 14 out of null spaces. 15 16 Logically Collective on MatNullSpace 17 18 Input Parameters: 19 + sp - the null space object 20 . rem - the function that removes the null space 21 - ctx - context for the remove function 22 23 Level: advanced 24 25 .keywords: PC, null space, create 26 27 .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceCreate() 28 @*/ 29 PetscErrorCode MatNullSpaceSetFunction(MatNullSpace sp, PetscErrorCode (*rem)(MatNullSpace,Vec,void*),void *ctx) 30 { 31 PetscFunctionBegin; 32 PetscValidHeaderSpecific(sp,MAT_NULLSPACE_CLASSID,1); 33 sp->remove = rem; 34 sp->rmctx = ctx; 35 PetscFunctionReturn(0); 36 } 37 38 #undef __FUNCT__ 39 #define __FUNCT__ "MatNullSpaceGetVecs" 40 /*@C 41 MatNullSpaceGetVecs - get vectors defining the null space 42 43 Not Collective 44 45 Input Arguments: 46 . sp - null space object 47 48 Output Arguments: 49 + has_cnst - PETSC_TRUE if the null space contains the constant vector, otherwise PETSC_FALSE 50 . n - number of vectors (excluding constant vector) in null space 51 - vecs - orthonormal vectors that span the null space (excluding the constant vector) 52 53 Level: developer 54 55 Notes: 56 These vectors and the array are owned by the MatNullSpace and should not be destroyed or freeded by the caller 57 58 .seealso: MatNullSpaceCreate(), MatGetNullSpace(), MatGetNearNullSpace() 59 @*/ 60 PetscErrorCode MatNullSpaceGetVecs(MatNullSpace sp,PetscBool *has_const,PetscInt *n,const Vec **vecs) 61 { 62 63 PetscFunctionBegin; 64 PetscValidHeaderSpecific(sp,MAT_NULLSPACE_CLASSID,1); 65 if (has_const) *has_const = sp->has_cnst; 66 if (n) *n = sp->n; 67 if (vecs) *vecs = sp->vecs; 68 PetscFunctionReturn(0); 69 } 70 71 #undef __FUNCT__ 72 #define __FUNCT__ "MatNullSpaceCreateRigidBody" 73 /*@ 74 MatNullSpaceCreateRigidBody - create rigid body modes from coordinates 75 76 Collective on Vec 77 78 Input Argument: 79 . coords - block of coordinates of each node, must have block size set 80 81 Output Argument: 82 . sp - the null space 83 84 Level: advanced 85 86 .seealso: MatNullSpaceCreate() 87 @*/ 88 PetscErrorCode MatNullSpaceCreateRigidBody(Vec coords,MatNullSpace *sp) 89 { 90 PetscErrorCode ierr; 91 const PetscScalar *x; 92 PetscScalar *v[6],dots[5]; 93 Vec vec[6]; 94 PetscInt n,N,dim,nmodes,i,j; 95 96 PetscFunctionBegin; 97 ierr = VecGetBlockSize(coords,&dim);CHKERRQ(ierr); 98 ierr = VecGetLocalSize(coords,&n);CHKERRQ(ierr); 99 ierr = VecGetSize(coords,&N);CHKERRQ(ierr); 100 n /= dim; 101 N /= dim; 102 switch (dim) { 103 case 1: 104 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_TRUE,0,NULL,sp);CHKERRQ(ierr); 105 break; 106 case 2: 107 case 3: 108 nmodes = (dim == 2) ? 3 : 6; 109 ierr = VecCreate(PetscObjectComm((PetscObject)coords),&vec[0]);CHKERRQ(ierr); 110 ierr = VecSetSizes(vec[0],dim*n,dim*N);CHKERRQ(ierr); 111 ierr = VecSetBlockSize(vec[0],dim);CHKERRQ(ierr); 112 ierr = VecSetUp(vec[0]);CHKERRQ(ierr); 113 for (i=1; i<nmodes; i++) {ierr = VecDuplicate(vec[0],&vec[i]);CHKERRQ(ierr);} 114 for (i=0; i<nmodes; i++) {ierr = VecGetArray(vec[i],&v[i]);CHKERRQ(ierr);} 115 ierr = VecGetArrayRead(coords,&x);CHKERRQ(ierr); 116 for (i=0; i<n; i++) { 117 if (dim == 2) { 118 v[0][i*2+0] = 1./N; 119 v[0][i*2+1] = 0.; 120 v[1][i*2+0] = 0.; 121 v[1][i*2+1] = 1./N; 122 /* Rotations */ 123 v[2][i*2+0] = -x[i*2+1]; 124 v[2][i*2+1] = x[i*2+0]; 125 } else { 126 v[0][i*3+0] = 1./N; 127 v[0][i*3+1] = 0.; 128 v[0][i*3+2] = 0.; 129 v[1][i*3+0] = 0.; 130 v[1][i*3+1] = 1./N; 131 v[1][i*3+2] = 0.; 132 v[2][i*3+0] = 0.; 133 v[2][i*3+1] = 0.; 134 v[2][i*3+2] = 1./N; 135 136 v[3][i*3+0] = x[i*3+1]; 137 v[3][i*3+1] = -x[i*3+0]; 138 v[3][i*3+2] = 0.; 139 v[4][i*3+0] = 0.; 140 v[4][i*3+1] = -x[i*3+2]; 141 v[4][i*3+2] = x[i*3+1]; 142 v[5][i*3+0] = x[i*3+2]; 143 v[5][i*3+1] = 0.; 144 v[5][i*3+2] = -x[i*3+0]; 145 } 146 } 147 for (i=0; i<nmodes; i++) {ierr = VecRestoreArray(vec[i],&v[i]);CHKERRQ(ierr);} 148 ierr = VecRestoreArrayRead(coords,&x);CHKERRQ(ierr); 149 for (i=dim; i<nmodes; i++) { 150 /* Orthonormalize vec[i] against vec[0:i-1] */ 151 ierr = VecMDot(vec[i],i,vec,dots);CHKERRQ(ierr); 152 for (j=0; j<i; j++) dots[j] *= -1.; 153 ierr = VecMAXPY(vec[i],i,dots,vec);CHKERRQ(ierr); 154 ierr = VecNormalize(vec[i],NULL);CHKERRQ(ierr); 155 } 156 ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)coords),PETSC_FALSE,nmodes,vec,sp);CHKERRQ(ierr); 157 for (i=0; i<nmodes; i++) {ierr = VecDestroy(&vec[i]);CHKERRQ(ierr);} 158 } 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "MatNullSpaceView" 164 /*@C 165 MatNullSpaceView - Visualizes a null space object. 166 167 Collective on MatNullSpace 168 169 Input Parameters: 170 + matnull - the null space 171 - viewer - visualization context 172 173 Level: advanced 174 175 Fortran Note: 176 This routine is not supported in Fortran. 177 178 .seealso: MatNullSpaceCreate(), PetscViewerASCIIOpen() 179 @*/ 180 PetscErrorCode MatNullSpaceView(MatNullSpace sp,PetscViewer viewer) 181 { 182 PetscErrorCode ierr; 183 PetscBool iascii; 184 185 PetscFunctionBegin; 186 PetscValidHeaderSpecific(sp,MAT_NULLSPACE_CLASSID,1); 187 if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)sp)); 188 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 189 PetscCheckSameComm(sp,1,viewer,2); 190 191 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 192 if (iascii) { 193 PetscViewerFormat format; 194 PetscInt i; 195 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 196 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)sp,viewer);CHKERRQ(ierr); 197 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 198 ierr = PetscViewerASCIIPrintf(viewer,"Contains %D vector%s%s\n",sp->n,sp->n==1 ? "" : "s",sp->has_cnst ? " and the constant" : "");CHKERRQ(ierr); 199 if (sp->remove) {ierr = PetscViewerASCIIPrintf(viewer,"Has user-provided removal function\n");CHKERRQ(ierr);} 200 if (!(format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL)) { 201 for (i=0; i<sp->n; i++) { 202 ierr = VecView(sp->vecs[i],viewer);CHKERRQ(ierr); 203 } 204 } 205 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 206 } 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "MatNullSpaceCreate" 212 /*@ 213 MatNullSpaceCreate - Creates a data structure used to project vectors 214 out of null spaces. 215 216 Collective on MPI_Comm 217 218 Input Parameters: 219 + comm - the MPI communicator associated with the object 220 . has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE 221 . n - number of vectors (excluding constant vector) in null space 222 - vecs - the vectors that span the null space (excluding the constant vector); 223 these vectors must be orthonormal. These vectors are NOT copied, so do not change them 224 after this call. You should free the array that you pass in and destroy the vectors (this will reduce the reference count 225 for them by one). 226 227 Output Parameter: 228 . SP - the null space context 229 230 Level: advanced 231 232 Notes: See MatNullSpaceSetFunction() as an alternative way of providing the null space information instead of setting vecs. 233 234 If has_cnst is PETSC_TRUE you do not need to pass a constant vector in as a fourth argument to this routine, nor do you 235 need to pass in a function that eliminates the constant function into MatNullSpaceSetFunction(). 236 237 Users manual sections: 238 . sec_singular 239 240 .keywords: PC, null space, create 241 242 .seealso: MatNullSpaceDestroy(), MatNullSpaceRemove(), MatSetNullSpace(), MatNullSpace, MatNullSpaceSetFunction() 243 @*/ 244 PetscErrorCode MatNullSpaceCreate(MPI_Comm comm,PetscBool has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP) 245 { 246 MatNullSpace sp; 247 PetscErrorCode ierr; 248 PetscInt i; 249 250 PetscFunctionBegin; 251 if (n < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of vectors (given %D) cannot be negative",n); 252 if (n) PetscValidPointer(vecs,4); 253 for (i=0; i<n; i++) PetscValidHeaderSpecific(vecs[i],VEC_CLASSID,4); 254 PetscValidPointer(SP,5); 255 256 *SP = NULL; 257 ierr = MatInitializePackage();CHKERRQ(ierr); 258 259 ierr = PetscHeaderCreate(sp,MAT_NULLSPACE_CLASSID,"MatNullSpace","Null space","Mat",comm,MatNullSpaceDestroy,MatNullSpaceView);CHKERRQ(ierr); 260 261 sp->has_cnst = has_cnst; 262 sp->n = n; 263 sp->vecs = 0; 264 sp->alpha = 0; 265 sp->remove = 0; 266 sp->rmctx = 0; 267 268 if (n) { 269 ierr = PetscMalloc1(n,&sp->vecs);CHKERRQ(ierr); 270 ierr = PetscMalloc1(n,&sp->alpha);CHKERRQ(ierr); 271 ierr = PetscLogObjectMemory((PetscObject)sp,n*(sizeof(Vec)+sizeof(PetscScalar)));CHKERRQ(ierr); 272 for (i=0; i<n; i++) { 273 ierr = PetscObjectReference((PetscObject)vecs[i]);CHKERRQ(ierr); 274 sp->vecs[i] = vecs[i]; 275 } 276 } 277 278 *SP = sp; 279 PetscFunctionReturn(0); 280 } 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "MatNullSpaceDestroy" 284 /*@ 285 MatNullSpaceDestroy - Destroys a data structure used to project vectors 286 out of null spaces. 287 288 Collective on MatNullSpace 289 290 Input Parameter: 291 . sp - the null space context to be destroyed 292 293 Level: advanced 294 295 .keywords: PC, null space, destroy 296 297 .seealso: MatNullSpaceCreate(), MatNullSpaceRemove(), MatNullSpaceSetFunction() 298 @*/ 299 PetscErrorCode MatNullSpaceDestroy(MatNullSpace *sp) 300 { 301 PetscErrorCode ierr; 302 303 PetscFunctionBegin; 304 if (!*sp) PetscFunctionReturn(0); 305 PetscValidHeaderSpecific((*sp),MAT_NULLSPACE_CLASSID,1); 306 if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; PetscFunctionReturn(0);} 307 308 ierr = VecDestroyVecs((*sp)->n,&(*sp)->vecs);CHKERRQ(ierr); 309 ierr = PetscFree((*sp)->alpha);CHKERRQ(ierr); 310 ierr = PetscHeaderDestroy(sp);CHKERRQ(ierr); 311 PetscFunctionReturn(0); 312 } 313 314 #undef __FUNCT__ 315 #define __FUNCT__ "MatNullSpaceRemove" 316 /*@C 317 MatNullSpaceRemove - Removes all the components of a null space from a vector. 318 319 Collective on MatNullSpace 320 321 Input Parameters: 322 + sp - the null space context 323 - vec - the vector from which the null space is to be removed 324 325 Level: advanced 326 327 .keywords: PC, null space, remove 328 329 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction() 330 @*/ 331 PetscErrorCode MatNullSpaceRemove(MatNullSpace sp,Vec vec) 332 { 333 PetscScalar sum; 334 PetscInt i,N; 335 PetscErrorCode ierr; 336 337 PetscFunctionBegin; 338 PetscValidHeaderSpecific(sp,MAT_NULLSPACE_CLASSID,1); 339 PetscValidHeaderSpecific(vec,VEC_CLASSID,2); 340 341 if (sp->has_cnst) { 342 ierr = VecGetSize(vec,&N);CHKERRQ(ierr); 343 if (N > 0) { 344 ierr = VecSum(vec,&sum);CHKERRQ(ierr); 345 sum = sum/((PetscScalar)(-1.0*N)); 346 ierr = VecShift(vec,sum);CHKERRQ(ierr); 347 } 348 } 349 350 if (sp->n) { 351 ierr = VecMDot(vec,sp->n,sp->vecs,sp->alpha);CHKERRQ(ierr); 352 for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i]; 353 ierr = VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);CHKERRQ(ierr); 354 } 355 356 if (sp->remove) { 357 ierr = (*sp->remove)(sp,vec,sp->rmctx);CHKERRQ(ierr); 358 } 359 PetscFunctionReturn(0); 360 } 361 362 #undef __FUNCT__ 363 #define __FUNCT__ "MatNullSpaceTest" 364 /*@ 365 MatNullSpaceTest - Tests if the claimed null space is really a 366 null space of a matrix 367 368 Collective on MatNullSpace 369 370 Input Parameters: 371 + sp - the null space context 372 - mat - the matrix 373 374 Output Parameters: 375 . isNull - PETSC_TRUE if the nullspace is valid for this matrix 376 377 Level: advanced 378 379 .keywords: PC, null space, remove 380 381 .seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction() 382 @*/ 383 PetscErrorCode MatNullSpaceTest(MatNullSpace sp,Mat mat,PetscBool *isNull) 384 { 385 PetscScalar sum; 386 PetscReal nrm; 387 PetscInt j,n,N; 388 PetscErrorCode ierr; 389 Vec l,r; 390 PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,consistent = PETSC_TRUE; 391 PetscViewer viewer; 392 393 PetscFunctionBegin; 394 PetscValidHeaderSpecific(sp,MAT_NULLSPACE_CLASSID,1); 395 PetscValidHeaderSpecific(mat,MAT_CLASSID,2); 396 n = sp->n; 397 ierr = PetscOptionsGetBool(NULL,"-mat_null_space_test_view",&flg1,NULL);CHKERRQ(ierr); 398 ierr = PetscOptionsGetBool(NULL,"-mat_null_space_test_view_draw",&flg2,NULL);CHKERRQ(ierr); 399 400 if (n) { 401 ierr = VecDuplicate(sp->vecs[0],&l);CHKERRQ(ierr); 402 } else { 403 ierr = MatCreateVecs(mat,&l,NULL);CHKERRQ(ierr); 404 } 405 406 ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)sp),&viewer);CHKERRQ(ierr); 407 if (sp->has_cnst) { 408 ierr = VecDuplicate(l,&r);CHKERRQ(ierr); 409 ierr = VecGetSize(l,&N);CHKERRQ(ierr); 410 sum = 1.0/N; 411 ierr = VecSet(l,sum);CHKERRQ(ierr); 412 ierr = MatMult(mat,l,r);CHKERRQ(ierr); 413 ierr = VecNorm(r,NORM_2,&nrm);CHKERRQ(ierr); 414 if (nrm >= 1.e-7) consistent = PETSC_FALSE; 415 if (flg1) { 416 if (consistent) { 417 ierr = PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are likely null vector");CHKERRQ(ierr); 418 } else { 419 ierr = PetscPrintf(PetscObjectComm((PetscObject)sp),"Constants are unlikely null vector ");CHKERRQ(ierr); 420 } 421 ierr = PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * 1/N || = %g\n",(double)nrm);CHKERRQ(ierr); 422 } 423 if (!consistent && flg1) {ierr = VecView(r,viewer);CHKERRQ(ierr);} 424 if (!consistent && flg2) {ierr = VecView(r,viewer);CHKERRQ(ierr);} 425 ierr = VecDestroy(&r);CHKERRQ(ierr); 426 } 427 428 for (j=0; j<n; j++) { 429 ierr = (*mat->ops->mult)(mat,sp->vecs[j],l);CHKERRQ(ierr); 430 ierr = VecNorm(l,NORM_2,&nrm);CHKERRQ(ierr); 431 if (nrm >= 1.e-7) consistent = PETSC_FALSE; 432 if (flg1) { 433 if (consistent) { 434 ierr = PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D is likely null vector",j);CHKERRQ(ierr); 435 } else { 436 ierr = PetscPrintf(PetscObjectComm((PetscObject)sp),"Null vector %D unlikely null vector ",j);CHKERRQ(ierr); 437 consistent = PETSC_FALSE; 438 } 439 ierr = PetscPrintf(PetscObjectComm((PetscObject)sp),"|| A * v[%D] || = %g\n",j,(double)nrm);CHKERRQ(ierr); 440 } 441 if (!consistent && flg1) {ierr = VecView(l,viewer);CHKERRQ(ierr);} 442 if (!consistent && flg2) {ierr = VecView(l,viewer);CHKERRQ(ierr);} 443 } 444 445 if (sp->remove) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot test a null space provided as a function with MatNullSpaceSetFunction()"); 446 ierr = VecDestroy(&l);CHKERRQ(ierr); 447 if (isNull) *isNull = consistent; 448 PetscFunctionReturn(0); 449 } 450 451