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