1 #include <tao.h> /*I "tao.h" I*/ 2 #include <petsc-private/matimpl.h> 3 #include <../src/tao/matrix/submatfree.h> 4 #include "tao_util.h" /*I "tao_util.h" I*/ 5 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "VecWhichEqual" 9 /*@ 10 VecWhichEqual - Creates an index set containing the indices 11 where the vectors Vec1 and Vec2 have identical elements. 12 13 Collective on Vec 14 15 Input Parameters: 16 . Vec1, Vec2 - the two vectors to compare 17 18 OutputParameter: 19 . S - The index set containing the indices i where vec1[i] == vec2[i] 20 21 Level: advanced 22 @*/ 23 PetscErrorCode VecWhichEqual(Vec Vec1, Vec Vec2, IS * S) 24 { 25 PetscErrorCode ierr; 26 PetscInt i,n_same = 0; 27 PetscInt n,low,high,low2,high2; 28 PetscInt *same = NULL; 29 PetscReal *v1,*v2; 30 MPI_Comm comm; 31 32 PetscFunctionBegin; 33 PetscValidHeaderSpecific(Vec1,VEC_CLASSID,1); 34 PetscValidHeaderSpecific(Vec2,VEC_CLASSID,2); 35 PetscCheckSameComm(Vec1,1,Vec2,2); 36 37 ierr = VecGetOwnershipRange(Vec1, &low, &high);CHKERRQ(ierr); 38 ierr = VecGetOwnershipRange(Vec2, &low2, &high2);CHKERRQ(ierr); 39 if ( low != low2 || high != high2 ) SETERRQ(PETSC_COMM_SELF,1,"Vectors must have identical layout"); 40 41 ierr = VecGetLocalSize(Vec1,&n);CHKERRQ(ierr); 42 if (n>0){ 43 if (Vec1 == Vec2){ 44 ierr = VecGetArray(Vec1,&v1);CHKERRQ(ierr); 45 v2=v1; 46 } else { 47 ierr = VecGetArray(Vec1,&v1);CHKERRQ(ierr); 48 ierr = VecGetArray(Vec2,&v2);CHKERRQ(ierr); 49 } 50 51 ierr = PetscMalloc1( n,&same );CHKERRQ(ierr); 52 53 for (i=0; i<n; i++){ 54 if (v1[i] == v2[i]) {same[n_same]=low+i; n_same++;} 55 } 56 57 if (Vec1 == Vec2){ 58 ierr = VecRestoreArray(Vec1,&v1);CHKERRQ(ierr); 59 } else { 60 ierr = VecRestoreArray(Vec1,&v1);CHKERRQ(ierr); 61 ierr = VecRestoreArray(Vec2,&v2);CHKERRQ(ierr); 62 } 63 } 64 ierr = PetscObjectGetComm((PetscObject)Vec1,&comm);CHKERRQ(ierr); 65 ierr = ISCreateGeneral(comm,n_same,same,PETSC_OWN_POINTER,S);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "VecWhichLessThan" 71 /*@ 72 VecWhichLessThan - Creates an index set containing the indices 73 where the vectors Vec1 < Vec2 74 75 Collective on S 76 77 Input Parameters: 78 . Vec1, Vec2 - the two vectors to compare 79 80 OutputParameter: 81 . S - The index set containing the indices i where vec1[i] < vec2[i] 82 83 Level: advanced 84 @*/ 85 PetscErrorCode VecWhichLessThan(Vec Vec1, Vec Vec2, IS * S) 86 { 87 PetscErrorCode ierr; 88 PetscInt i; 89 PetscInt n,low,high,low2,high2,n_lt=0; 90 PetscInt *lt = NULL; 91 PetscReal *v1,*v2; 92 MPI_Comm comm; 93 94 PetscFunctionBegin; 95 PetscValidHeaderSpecific(Vec1,VEC_CLASSID,1); 96 PetscValidHeaderSpecific(Vec2,VEC_CLASSID,2); 97 PetscCheckSameComm(Vec1,1,Vec2,2); 98 99 ierr = VecGetOwnershipRange(Vec1, &low, &high);CHKERRQ(ierr); 100 ierr = VecGetOwnershipRange(Vec2, &low2, &high2);CHKERRQ(ierr); 101 if ( low != low2 || high != high2 ) SETERRQ(PETSC_COMM_SELF,1,"Vectors must haveidentical layout"); 102 103 ierr = VecGetLocalSize(Vec1,&n);CHKERRQ(ierr); 104 if (n>0){ 105 if (Vec1 == Vec2){ 106 ierr = VecGetArray(Vec1,&v1);CHKERRQ(ierr); 107 v2=v1; 108 } else { 109 ierr = VecGetArray(Vec1,&v1);CHKERRQ(ierr); 110 ierr = VecGetArray(Vec2,&v2);CHKERRQ(ierr); 111 } 112 ierr = PetscMalloc1(n,< );CHKERRQ(ierr); 113 114 for (i=0; i<n; i++){ 115 if (v1[i] < v2[i]) {lt[n_lt]=low+i; n_lt++;} 116 } 117 118 if (Vec1 == Vec2){ 119 ierr = VecRestoreArray(Vec1,&v1);CHKERRQ(ierr); 120 } else { 121 ierr = VecRestoreArray(Vec1,&v1);CHKERRQ(ierr); 122 ierr = VecRestoreArray(Vec2,&v2);CHKERRQ(ierr); 123 } 124 } 125 ierr = PetscObjectGetComm((PetscObject)Vec1,&comm);CHKERRQ(ierr); 126 ierr = ISCreateGeneral(comm,n_lt,lt,PETSC_OWN_POINTER,S);CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "VecWhichGreaterThan" 132 /*@ 133 VecWhichGreaterThan - Creates an index set containing the indices 134 where the vectors Vec1 > Vec2 135 136 Collective on S 137 138 Input Parameters: 139 . Vec1, Vec2 - the two vectors to compare 140 141 OutputParameter: 142 . S - The index set containing the indices i where vec1[i] > vec2[i] 143 144 Level: advanced 145 @*/ 146 PetscErrorCode VecWhichGreaterThan(Vec Vec1, Vec Vec2, IS * S) 147 { 148 PetscErrorCode ierr; 149 PetscInt n,low,high,low2,high2,n_gt=0,i; 150 PetscInt *gt=NULL; 151 PetscReal *v1,*v2; 152 MPI_Comm comm; 153 154 PetscFunctionBegin; 155 PetscValidHeaderSpecific(Vec1,VEC_CLASSID,1); 156 PetscValidHeaderSpecific(Vec2,VEC_CLASSID,2); 157 PetscCheckSameComm(Vec1,1,Vec2,2); 158 159 ierr = VecGetOwnershipRange(Vec1, &low, &high);CHKERRQ(ierr); 160 ierr = VecGetOwnershipRange(Vec2, &low2, &high2);CHKERRQ(ierr); 161 if ( low != low2 || high != high2 ) SETERRQ(PETSC_COMM_SELF,1,"Vectors must be have identical layout"); 162 163 ierr = VecGetLocalSize(Vec1,&n);CHKERRQ(ierr); 164 165 if (n>0){ 166 167 if (Vec1 == Vec2){ 168 ierr = VecGetArray(Vec1,&v1);CHKERRQ(ierr); 169 v2=v1; 170 } else { 171 ierr = VecGetArray(Vec1,&v1);CHKERRQ(ierr); 172 ierr = VecGetArray(Vec2,&v2);CHKERRQ(ierr); 173 } 174 175 ierr = PetscMalloc1(n, > );CHKERRQ(ierr); 176 177 for (i=0; i<n; i++){ 178 if (v1[i] > v2[i]) {gt[n_gt]=low+i; n_gt++;} 179 } 180 181 if (Vec1 == Vec2){ 182 ierr = VecRestoreArray(Vec1,&v1);CHKERRQ(ierr); 183 } else { 184 ierr = VecRestoreArray(Vec1,&v1);CHKERRQ(ierr); 185 ierr = VecRestoreArray(Vec2,&v2);CHKERRQ(ierr); 186 } 187 } 188 ierr = PetscObjectGetComm((PetscObject)Vec1,&comm);CHKERRQ(ierr); 189 ierr = ISCreateGeneral(comm,n_gt,gt,PETSC_OWN_POINTER,S);CHKERRQ(ierr); 190 PetscFunctionReturn(0); 191 } 192 193 #undef __FUNCT__ 194 #define __FUNCT__ "VecWhichBetween" 195 /*@ 196 VecWhichBetween - Creates an index set containing the indices 197 where VecLow < V < VecHigh 198 199 Collective on S 200 201 Input Parameters: 202 + VecLow - lower bound 203 . V - Vector to compare 204 - VecHigh - higher bound 205 206 OutputParameter: 207 . S - The index set containing the indices i where veclow[i] < v[i] < vechigh[i] 208 209 Level: advanced 210 @*/ 211 PetscErrorCode VecWhichBetween(Vec VecLow, Vec V, Vec VecHigh, IS *S) 212 { 213 214 PetscErrorCode ierr; 215 PetscInt n,low,high,low2,high2,low3,high3,n_vm=0; 216 PetscInt *vm,i; 217 PetscReal *v1,*v2,*vmiddle; 218 MPI_Comm comm; 219 220 PetscValidHeaderSpecific(V,VEC_CLASSID,2); 221 PetscCheckSameComm(V,2,VecLow,1); PetscCheckSameComm(V,2,VecHigh,3); 222 223 ierr = VecGetOwnershipRange(VecLow, &low, &high);CHKERRQ(ierr); 224 ierr = VecGetOwnershipRange(VecHigh, &low2, &high2);CHKERRQ(ierr); 225 ierr = VecGetOwnershipRange(V, &low3, &high3);CHKERRQ(ierr); 226 if ( low!=low2 || high!=high2 || low!=low3 || high!=high3) SETERRQ(PETSC_COMM_SELF,1,"Vectors must have identical layout"); 227 228 ierr = VecGetLocalSize(VecLow,&n);CHKERRQ(ierr); 229 if (n>0){ 230 ierr = VecGetArray(VecLow,&v1);CHKERRQ(ierr); 231 if (VecLow != VecHigh){ 232 ierr = VecGetArray(VecHigh,&v2);CHKERRQ(ierr); 233 } else { 234 v2=v1; 235 } 236 if ( V != VecLow && V != VecHigh){ 237 ierr = VecGetArray(V,&vmiddle);CHKERRQ(ierr); 238 } else if ( V==VecLow ){ 239 vmiddle=v1; 240 } else { 241 vmiddle =v2; 242 } 243 244 ierr = PetscMalloc1(n, &vm );CHKERRQ(ierr); 245 246 for (i=0; i<n; i++){ 247 if (v1[i] < vmiddle[i] && vmiddle[i] < v2[i]) {vm[n_vm]=low+i; n_vm++;} 248 } 249 250 ierr = VecRestoreArray(VecLow,&v1);CHKERRQ(ierr); 251 if (VecLow != VecHigh){ 252 ierr = VecRestoreArray(VecHigh,&v2);CHKERRQ(ierr); 253 } 254 if ( V != VecLow && V != VecHigh){ 255 ierr = VecRestoreArray(V,&vmiddle);CHKERRQ(ierr); 256 } 257 } 258 ierr = PetscObjectGetComm((PetscObject)V,&comm);CHKERRQ(ierr); 259 ierr = ISCreateGeneral(comm,n_vm,vm,PETSC_OWN_POINTER,S);CHKERRQ(ierr); 260 PetscFunctionReturn(0); 261 } 262 263 264 #undef __FUNCT__ 265 #define __FUNCT__ "VecWhichBetweenOrEqual" 266 /*@ 267 VecWhichBetweenOrEqual - Creates an index set containing the indices 268 where VecLow <= V <= VecHigh 269 270 Collective on S 271 272 Input Parameters: 273 + VecLow - lower bound 274 . V - Vector to compare 275 - VecHigh - higher bound 276 277 OutputParameter: 278 . S - The index set containing the indices i where veclow[i] <= v[i] <= vechigh[i] 279 280 Level: advanced 281 @*/ 282 283 PetscErrorCode VecWhichBetweenOrEqual(Vec VecLow, Vec V, Vec VecHigh, IS * S) 284 { 285 PetscErrorCode ierr; 286 PetscInt n,low,high,low2,high2,low3,high3,n_vm=0,i; 287 PetscInt *vm = NULL; 288 PetscReal *v1,*v2,*vmiddle; 289 MPI_Comm comm; 290 291 PetscValidHeaderSpecific(V,VEC_CLASSID,2); 292 PetscCheckSameComm(V,2,VecLow,1); PetscCheckSameComm(V,2,VecHigh,3); 293 294 ierr = VecGetOwnershipRange(VecLow, &low, &high);CHKERRQ(ierr); 295 ierr = VecGetOwnershipRange(VecHigh, &low2, &high2);CHKERRQ(ierr); 296 ierr = VecGetOwnershipRange(V, &low3, &high3);CHKERRQ(ierr); 297 if ( low!=low2 || high!=high2 || low!=low3 || high!=high3 ) SETERRQ(PETSC_COMM_SELF,1,"Vectors must have identical layout"); 298 299 ierr = VecGetLocalSize(VecLow,&n);CHKERRQ(ierr); 300 301 if (n>0){ 302 ierr = VecGetArray(VecLow,&v1);CHKERRQ(ierr); 303 if (VecLow != VecHigh){ 304 ierr = VecGetArray(VecHigh,&v2);CHKERRQ(ierr); 305 } else { 306 v2=v1; 307 } 308 if ( V != VecLow && V != VecHigh){ 309 ierr = VecGetArray(V,&vmiddle);CHKERRQ(ierr); 310 } else if ( V==VecLow ){ 311 vmiddle=v1; 312 } else { 313 vmiddle =v2; 314 } 315 316 ierr = PetscMalloc1(n, &vm );CHKERRQ(ierr); 317 318 for (i=0; i<n; i++){ 319 if (v1[i] <= vmiddle[i] && vmiddle[i] <= v2[i]) {vm[n_vm]=low+i; n_vm++;} 320 } 321 322 ierr = VecRestoreArray(VecLow,&v1);CHKERRQ(ierr); 323 if (VecLow != VecHigh){ 324 ierr = VecRestoreArray(VecHigh,&v2);CHKERRQ(ierr); 325 } 326 if ( V != VecLow && V != VecHigh){ 327 ierr = VecRestoreArray(V,&vmiddle);CHKERRQ(ierr); 328 } 329 } 330 ierr = PetscObjectGetComm((PetscObject)V,&comm);CHKERRQ(ierr); 331 ierr = ISCreateGeneral(comm,n_vm,vm,PETSC_OWN_POINTER,S);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 336 #undef __FUNCT__ 337 #define __FUNCT__ "VecGetSubVec" 338 /*@ 339 VecGetSubVec - Gets a subvector using the IS 340 341 Input Parameters: 342 + vfull - the full matrix 343 . is - the index set for the subvector 344 . reduced_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK, TAO_SUBSET_MATRIXFREE) 345 - maskvalue - the value to set the unused vector elements to (for TAO_SUBSET_MASK or TAO_SUBSET_MATRIXFREE) 346 347 Output Parameters: 348 . vreduced - the subvector 349 350 Note: 351 maskvalue should usually be 0.0, unless a pointwise divide will be used. 352 @*/ 353 PetscErrorCode VecGetSubVec(Vec vfull, IS is, PetscInt reduced_type, PetscReal maskvalue, Vec *vreduced) 354 { 355 PetscErrorCode ierr; 356 PetscInt nfull,nreduced,nreduced_local,rlow,rhigh,flow,fhigh; 357 PetscInt i,nlocal; 358 PetscReal *fv,*rv; 359 const PetscInt *s; 360 IS ident; 361 VecType vtype; 362 VecScatter scatter; 363 MPI_Comm comm; 364 365 PetscFunctionBegin; 366 PetscValidHeaderSpecific(vfull,VEC_CLASSID,1); 367 PetscValidHeaderSpecific(is,IS_CLASSID,2); 368 369 ierr = VecGetSize(vfull, &nfull);CHKERRQ(ierr); 370 ierr = ISGetSize(is, &nreduced);CHKERRQ(ierr); 371 372 if (nreduced == nfull) { 373 ierr = VecDestroy(vreduced);CHKERRQ(ierr); 374 ierr = VecDuplicate(vfull,vreduced);CHKERRQ(ierr); 375 ierr = VecCopy(vfull,*vreduced);CHKERRQ(ierr); 376 } else { 377 switch (reduced_type) { 378 case TAO_SUBSET_SUBVEC: 379 ierr = VecGetType(vfull,&vtype);CHKERRQ(ierr); 380 ierr = VecGetOwnershipRange(vfull,&flow,&fhigh);CHKERRQ(ierr); 381 ierr = ISGetLocalSize(is,&nreduced_local);CHKERRQ(ierr); 382 ierr = PetscObjectGetComm((PetscObject)vfull,&comm);CHKERRQ(ierr); 383 if (*vreduced) { 384 ierr = VecDestroy(vreduced);CHKERRQ(ierr); 385 } 386 ierr = VecCreate(comm,vreduced);CHKERRQ(ierr); 387 ierr = VecSetType(*vreduced,vtype);CHKERRQ(ierr); 388 389 ierr = VecSetSizes(*vreduced,nreduced_local,nreduced);CHKERRQ(ierr); 390 ierr = VecGetOwnershipRange(*vreduced,&rlow,&rhigh);CHKERRQ(ierr); 391 ierr = ISCreateStride(comm,nreduced_local,rlow,1,&ident);CHKERRQ(ierr); 392 ierr = VecScatterCreate(vfull,is,*vreduced,ident,&scatter);CHKERRQ(ierr); 393 ierr = VecScatterBegin(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 394 ierr = VecScatterEnd(scatter,vfull,*vreduced,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 395 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 396 ierr = ISDestroy(&ident);CHKERRQ(ierr); 397 break; 398 399 case TAO_SUBSET_MASK: 400 case TAO_SUBSET_MATRIXFREE: 401 /* vr[i] = vf[i] if i in is 402 vr[i] = 0 otherwise */ 403 if (*vreduced == NULL) { 404 ierr = VecDuplicate(vfull,vreduced);CHKERRQ(ierr); 405 } 406 407 ierr = VecSet(*vreduced,maskvalue);CHKERRQ(ierr); 408 ierr = ISGetLocalSize(is,&nlocal);CHKERRQ(ierr); 409 ierr = VecGetOwnershipRange(vfull,&flow,&fhigh);CHKERRQ(ierr); 410 ierr = VecGetArray(vfull,&fv);CHKERRQ(ierr); 411 ierr = VecGetArray(*vreduced,&rv);CHKERRQ(ierr); 412 ierr = ISGetIndices(is,&s);CHKERRQ(ierr); 413 if (nlocal > (fhigh-flow)) SETERRQ2(PETSC_COMM_WORLD,1,"IS local size %d > Vec local size %d",nlocal,fhigh-flow); 414 for (i=0;i<nlocal;i++) { 415 rv[s[i]-flow] = fv[s[i]-flow]; 416 } 417 ierr = ISRestoreIndices(is,&s);CHKERRQ(ierr); 418 ierr = VecRestoreArray(vfull,&fv);CHKERRQ(ierr); 419 ierr = VecRestoreArray(*vreduced,&rv);CHKERRQ(ierr); 420 break; 421 } 422 } 423 PetscFunctionReturn(0); 424 } 425 426 #undef __FUNCT__ 427 #define __FUNCT__ "VecReducedXPY" 428 /*@ 429 VecReducedXPY - Adds a reduced vector to the appropriate elements of a full-space vector. 430 431 Input Parameters: 432 + vfull - the full-space vector 433 . vreduced - the reduced-space vector 434 - is - the index set for the reduced space 435 436 Output Parameters: 437 . vfull - the sum of the full-space vector and reduced-space vector 438 @*/ 439 PetscErrorCode VecReducedXPY(Vec vfull, Vec vreduced, IS is) 440 { 441 VecScatter scatter; 442 IS ident; 443 PetscInt nfull,nreduced,rlow,rhigh; 444 MPI_Comm comm; 445 PetscErrorCode ierr; 446 447 PetscFunctionBegin; 448 PetscValidHeaderSpecific(vfull,VEC_CLASSID,1); 449 PetscValidHeaderSpecific(vreduced,VEC_CLASSID,2); 450 PetscValidHeaderSpecific(is,IS_CLASSID,3); 451 ierr = VecGetSize(vfull,&nfull);CHKERRQ(ierr); 452 ierr = VecGetSize(vreduced,&nreduced);CHKERRQ(ierr); 453 454 if (nfull == nreduced) { /* Also takes care of masked vectors */ 455 ierr = VecAXPY(vfull,1.0,vreduced);CHKERRQ(ierr); 456 } else { 457 ierr = PetscObjectGetComm((PetscObject)vfull,&comm);CHKERRQ(ierr); 458 ierr = VecGetOwnershipRange(vreduced,&rlow,&rhigh);CHKERRQ(ierr); 459 ierr = ISCreateStride(comm,rhigh-rlow,rlow,1,&ident);CHKERRQ(ierr); 460 ierr = VecScatterCreate(vreduced,ident,vfull,is,&scatter);CHKERRQ(ierr); 461 ierr = VecScatterBegin(scatter,vreduced,vfull,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 462 ierr = VecScatterEnd(scatter,vreduced,vfull,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 463 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr); 464 ierr = ISDestroy(&ident);CHKERRQ(ierr); 465 } 466 PetscFunctionReturn(0); 467 } 468 469 470 #undef __FUNCT__ 471 #define __FUNCT__ "ISCreateComplement" 472 /*@ 473 ISCreateComplement - Creates the complement of the the index set 474 475 Collective on IS 476 477 Input Parameter: 478 + S - a PETSc IS 479 - V - the reference vector space 480 481 Output Parameter: 482 . T - the complement of S 483 484 485 .seealso ISCreateGeneral() 486 487 Level: advanced 488 @*/ 489 PetscErrorCode ISCreateComplement(IS S, Vec V, IS *T) 490 { 491 PetscErrorCode ierr; 492 PetscInt i,nis,nloc,high,low,n=0; 493 const PetscInt *s; 494 PetscInt *tt,*ss; 495 MPI_Comm comm; 496 497 PetscFunctionBegin; 498 PetscValidHeaderSpecific(S,IS_CLASSID,1); 499 PetscValidHeaderSpecific(V,VEC_CLASSID,2); 500 501 ierr = VecGetOwnershipRange(V,&low,&high);CHKERRQ(ierr); 502 ierr = VecGetLocalSize(V,&nloc);CHKERRQ(ierr); 503 ierr = ISGetLocalSize(S,&nis);CHKERRQ(ierr); 504 ierr = PetscMalloc1(nloc,&tt );CHKERRQ(ierr); 505 ierr = PetscMalloc1(nloc,&ss );CHKERRQ(ierr); 506 507 ierr = ISGetIndices(S, &s);CHKERRQ(ierr); 508 for (i=low; i<high; i++){ tt[i-low]=i; } 509 for (i=0; i<nis; i++){ tt[s[i]-low] = -2; } 510 511 for (i=0; i<nloc; i++){ 512 if (tt[i]>-1){ ss[n]=tt[i]; n++; } 513 } 514 ierr = ISRestoreIndices(S, &s);CHKERRQ(ierr); 515 516 ierr = PetscObjectGetComm((PetscObject)S,&comm);CHKERRQ(ierr); 517 ierr = ISCreateGeneral(comm,n,ss,PETSC_COPY_VALUES,T);CHKERRQ(ierr); 518 ierr = PetscFree(tt);CHKERRQ(ierr); 519 ierr = PetscFree(ss);CHKERRQ(ierr); 520 PetscFunctionReturn(0); 521 } 522 523 #undef __FUNCT__ 524 #define __FUNCT__ "VecISSetToConstant" 525 /*@ 526 VecISSetToConstant - Sets the elements of a vector, specified by an index set, to a constant 527 528 Input Parameter: 529 + S - a PETSc IS 530 . c - the constant 531 - V - a Vec 532 533 .seealso VecSet() 534 535 Level: advanced 536 @*/ 537 PetscErrorCode VecISSetToConstant(IS S, PetscReal c, Vec V) 538 { 539 PetscErrorCode ierr; 540 PetscInt nloc,low,high,i; 541 const PetscInt *s; 542 PetscReal *v; 543 544 PetscFunctionBegin; 545 PetscValidHeaderSpecific(V,VEC_CLASSID,3); 546 PetscValidHeaderSpecific(S,IS_CLASSID,1); 547 PetscValidType(V,3); 548 PetscCheckSameComm(V,3,S,1); 549 550 ierr = VecGetOwnershipRange(V, &low, &high);CHKERRQ(ierr); 551 ierr = ISGetLocalSize(S,&nloc);CHKERRQ(ierr); 552 ierr = ISGetIndices(S, &s);CHKERRQ(ierr); 553 ierr = VecGetArray(V,&v);CHKERRQ(ierr); 554 for (i=0; i<nloc; i++){ 555 v[s[i]-low] = c; 556 } 557 ierr = ISRestoreIndices(S, &s);CHKERRQ(ierr); 558 ierr = VecRestoreArray(V,&v);CHKERRQ(ierr); 559 PetscFunctionReturn(0); 560 } 561 562 #undef __FUNCT__ 563 #define __FUNCT__ "MatGetSubMat" 564 /*@ 565 MatGetSubMat - Gets a submatrix using the IS 566 567 Input Parameters: 568 + M - the full matrix (n x n) 569 . is - the index set for the submatrix (both row and column index sets need to be the same) 570 . v1 - work vector of dimension n, needed for TAO_SUBSET_MASK option 571 - subset_type - the method TAO is using for subsetting (TAO_SUBSET_SUBVEC, TAO_SUBSET_MASK, 572 TAO_SUBSET_MATRIXFREE) 573 574 Output Parameters: 575 . Msub - the submatrix 576 @*/ 577 PetscErrorCode MatGetSubMat(Mat M, IS is, Vec v1, TaoSubsetType subset_type, Mat *Msub) 578 { 579 PetscErrorCode ierr; 580 IS iscomp; 581 PetscBool flg; 582 583 PetscFunctionBegin; 584 PetscValidHeaderSpecific(M,MAT_CLASSID,1); 585 PetscValidHeaderSpecific(is,IS_CLASSID,2); 586 ierr = MatDestroy(Msub);CHKERRQ(ierr); 587 switch (subset_type) { 588 case TAO_SUBSET_SUBVEC: 589 ierr = MatGetSubMatrix(M, is, is, MAT_INITIAL_MATRIX, Msub);CHKERRQ(ierr); 590 break; 591 592 case TAO_SUBSET_MASK: 593 /* Get Reduced Hessian 594 Msub[i,j] = M[i,j] if i,j in Free_Local or i==j 595 Msub[i,j] = 0 if i!=j and i or j not in Free_Local 596 */ 597 ierr = PetscOptionsBool("-different_submatrix","use separate hessian matrix when computing submatrices","TaoSubsetType",PETSC_FALSE,&flg,NULL); 598 if (flg == PETSC_TRUE) { 599 ierr = MatDuplicate(M, MAT_COPY_VALUES, Msub);CHKERRQ(ierr); 600 } else { 601 /* Act on hessian directly (default) */ 602 ierr = PetscObjectReference((PetscObject)M);CHKERRQ(ierr); 603 *Msub = M; 604 } 605 /* Save the diagonal to temporary vector */ 606 ierr = MatGetDiagonal(*Msub,v1);CHKERRQ(ierr); 607 608 /* Zero out rows and columns */ 609 ierr = ISCreateComplement(is,v1,&iscomp);CHKERRQ(ierr); 610 611 /* Use v1 instead of 0 here because of PETSc bug */ 612 ierr = MatZeroRowsColumnsIS(*Msub,iscomp,1.0,v1,v1);CHKERRQ(ierr); 613 614 ierr = ISDestroy(&iscomp);CHKERRQ(ierr); 615 break; 616 case TAO_SUBSET_MATRIXFREE: 617 ierr = ISCreateComplement(is,v1,&iscomp);CHKERRQ(ierr); 618 ierr = MatCreateSubMatrixFree(M,iscomp,iscomp,Msub);CHKERRQ(ierr); 619 ierr = ISDestroy(&iscomp);CHKERRQ(ierr); 620 break; 621 } 622 PetscFunctionReturn(0); 623 } 624