1# -------------------------------------------------------------------- 2 3class VecType(object): 4 SEQ = S_(VECSEQ) 5 MPI = S_(VECMPI) 6 STANDARD = S_(VECSTANDARD) 7 SHARED = S_(VECSHARED) 8 SEQVIENNACL= S_(VECSEQVIENNACL) 9 MPIVIENNACL= S_(VECMPIVIENNACL) 10 VIENNACL = S_(VECVIENNACL) 11 SEQCUDA = S_(VECSEQCUDA) 12 MPICUDA = S_(VECMPICUDA) 13 CUDA = S_(VECCUDA) 14 NEST = S_(VECNEST) 15 NODE = S_(VECNODE) 16 SEQKOKKOS = S_(VECSEQKOKKOS) 17 MPIKOKKOS = S_(VECMPIKOKKOS) 18 KOKKOS = S_(VECKOKKOS) 19 20class VecOption(object): 21 IGNORE_OFF_PROC_ENTRIES = VEC_IGNORE_OFF_PROC_ENTRIES 22 IGNORE_NEGATIVE_INDICES = VEC_IGNORE_NEGATIVE_INDICES 23 24# -------------------------------------------------------------------- 25 26cdef class Vec(Object): 27 28 Type = VecType 29 Option = VecOption 30 31 # 32 33 def __cinit__(self): 34 self.obj = <PetscObject*> &self.vec 35 self.vec = NULL 36 37 # unary operations 38 39 def __pos__(self): 40 return vec_pos(self) 41 42 def __neg__(self): 43 return vec_neg(self) 44 45 def __abs__(self): 46 return vec_abs(self) 47 48 # inplace binary operations 49 50 def __iadd__(self, other): 51 return vec_iadd(self, other) 52 53 def __isub__(self, other): 54 return vec_isub(self, other) 55 56 def __imul__(self, other): 57 return vec_imul(self, other) 58 59 def __idiv__(self, other): 60 return vec_idiv(self, other) 61 62 def __itruediv__(self, other): 63 return vec_idiv(self, other) 64 65 # binary operations 66 67 def __add__(self, other): 68 if isinstance(self, Vec): 69 return vec_add(self, other) 70 else: 71 return vec_radd(other, self) 72 73 def __sub__(self, other): 74 if isinstance(self, Vec): 75 return vec_sub(self, other) 76 else: 77 return vec_rsub(other, self) 78 79 def __mul__(self, other): 80 if isinstance(self, Vec): 81 return vec_mul(self, other) 82 else: 83 return vec_rmul(other, self) 84 85 def __div__(self, other): 86 if isinstance(self, Vec): 87 return vec_div(self, other) 88 else: 89 return vec_rdiv(other, self) 90 91 def __truediv__(self, other): 92 if isinstance(self, Vec): 93 return vec_div(self, other) 94 else: 95 return vec_rdiv(other, self) 96 97 # 98 99 #def __len__(self): 100 # cdef PetscInt size = 0 101 # CHKERR( VecGetSize(self.vec, &size) ) 102 # return <Py_ssize_t>size 103 104 def __getitem__(self, i): 105 return vec_getitem(self, i) 106 107 def __setitem__(self, i, v): 108 vec_setitem(self, i, v) 109 110 # buffer interface (PEP 3118) 111 112 def __getbuffer__(self, Py_buffer *view, int flags): 113 cdef _Vec_buffer buf = _Vec_buffer(self) 114 buf.acquirebuffer(view, flags) 115 116 def __releasebuffer__(self, Py_buffer *view): 117 cdef _Vec_buffer buf = <_Vec_buffer>(view.obj) 118 buf.releasebuffer(view) 119 <void>self # unused 120 121 # 'with' statement (PEP 343) 122 123 def __enter__(self): 124 cdef _Vec_buffer buf = _Vec_buffer(self) 125 self.set_attr('__buffer__', buf) 126 return buf.enter() 127 128 def __exit__(self, *exc): 129 cdef _Vec_buffer buf = self.get_attr('__buffer__') 130 self.set_attr('__buffer__', None) 131 return buf.exit() 132 133 # 134 135 def view(self, Viewer viewer=None): 136 cdef PetscViewer vwr = NULL 137 if viewer is not None: vwr = viewer.vwr 138 CHKERR( VecView(self.vec, vwr) ) 139 140 def destroy(self): 141 CHKERR( VecDestroy(&self.vec) ) 142 return self 143 144 def create(self, comm=None): 145 cdef MPI_Comm ccomm = def_Comm(comm, PETSC_COMM_DEFAULT) 146 cdef PetscVec newvec = NULL 147 CHKERR( VecCreate(ccomm, &newvec) ) 148 PetscCLEAR(self.obj); self.vec = newvec 149 return self 150 151 def setType(self, vec_type): 152 cdef PetscVecType cval = NULL 153 vec_type = str2bytes(vec_type, &cval) 154 CHKERR( VecSetType(self.vec, cval) ) 155 156 def setSizes(self, size, bsize=None): 157 cdef PetscInt bs=0, n=0, N=0 158 Vec_Sizes(size, bsize, &bs, &n, &N) 159 CHKERR( VecSetSizes(self.vec, n, N) ) 160 if bs != PETSC_DECIDE: 161 CHKERR( VecSetBlockSize(self.vec, bs) ) 162 163 # 164 165 def createSeq(self, size, bsize=None, comm=None): 166 cdef MPI_Comm ccomm = def_Comm(comm, PETSC_COMM_SELF) 167 cdef PetscInt bs=0, n=0, N=0 168 Vec_Sizes(size, bsize, &bs, &n, &N) 169 Sys_Layout(ccomm, bs, &n, &N) 170 if bs == PETSC_DECIDE: bs = 1 171 cdef PetscVec newvec = NULL 172 CHKERR( VecCreate(ccomm,&newvec) ) 173 CHKERR( VecSetSizes(newvec, n, N) ) 174 CHKERR( VecSetBlockSize(newvec, bs) ) 175 CHKERR( VecSetType(newvec, VECSEQ) ) 176 PetscCLEAR(self.obj); self.vec = newvec 177 return self 178 179 def createMPI(self, size, bsize=None, comm=None): 180 cdef MPI_Comm ccomm = def_Comm(comm, PETSC_COMM_DEFAULT) 181 cdef PetscInt bs=0, n=0, N=0 182 Vec_Sizes(size, bsize, &bs, &n, &N) 183 Sys_Layout(ccomm, bs, &n, &N) 184 if bs == PETSC_DECIDE: bs = 1 185 cdef PetscVec newvec = NULL 186 CHKERR( VecCreate(ccomm, &newvec) ) 187 CHKERR( VecSetSizes(newvec, n, N) ) 188 CHKERR( VecSetBlockSize(newvec, bs) ) 189 CHKERR( VecSetType(newvec, VECMPI) ) 190 PetscCLEAR(self.obj); self.vec = newvec 191 return self 192 193 def createWithArray(self, array, size=None, bsize=None, comm=None): 194 cdef PetscInt na=0 195 cdef PetscScalar *sa=NULL 196 array = iarray_s(array, &na, &sa) 197 if size is None: size = (toInt(na), toInt(PETSC_DECIDE)) 198 cdef MPI_Comm ccomm = def_Comm(comm, PETSC_COMM_DEFAULT) 199 cdef PetscInt bs=0, n=0, N=0 200 Vec_Sizes(size, bsize, &bs, &n, &N) 201 Sys_Layout(ccomm, bs, &n, &N) 202 if bs == PETSC_DECIDE: bs = 1 203 if na < n: raise ValueError( 204 "array size %d and vector local size %d block size %d" % 205 (toInt(na), toInt(n), toInt(bs))) 206 cdef PetscVec newvec = NULL 207 if comm_size(ccomm) == 1: 208 CHKERR( VecCreateSeqWithArray(ccomm,bs,N,sa,&newvec) ) 209 else: 210 CHKERR( VecCreateMPIWithArray(ccomm,bs,n,N,sa,&newvec) ) 211 PetscCLEAR(self.obj); self.vec = newvec 212 self.set_attr('__array__', array) 213 return self 214 215 def createGhost(self, ghosts, size, bsize=None, comm=None): 216 cdef MPI_Comm ccomm = def_Comm(comm, PETSC_COMM_DEFAULT) 217 cdef PetscInt ng=0, *ig=NULL 218 ghosts = iarray_i(ghosts, &ng, &ig) 219 cdef PetscInt bs=0, n=0, N=0 220 Vec_Sizes(size, bsize, &bs, &n, &N) 221 Sys_Layout(ccomm, bs, &n, &N) 222 cdef PetscVec newvec = NULL 223 if bs == PETSC_DECIDE: 224 CHKERR( VecCreateGhost( 225 ccomm, n, N, ng, ig, &newvec) ) 226 else: 227 CHKERR( VecCreateGhostBlock( 228 ccomm, bs, n, N, ng, ig, &newvec) ) 229 PetscCLEAR(self.obj); self.vec = newvec 230 return self 231 232 def createGhostWithArray(self, ghosts, array, 233 size=None, bsize=None, comm=None): 234 cdef MPI_Comm ccomm = def_Comm(comm, PETSC_COMM_DEFAULT) 235 cdef PetscInt ng=0, *ig=NULL 236 ghosts = iarray_i(ghosts, &ng, &ig) 237 cdef PetscInt na=0 238 cdef PetscScalar *sa=NULL 239 array = oarray_s(array, &na, &sa) 240 cdef PetscInt b = 1 if bsize is None else asInt(bsize) 241 if size is None: size = (toInt(na-ng*b), toInt(PETSC_DECIDE)) 242 cdef PetscInt bs=0, n=0, N=0 243 Vec_Sizes(size, bsize, &bs, &n, &N) 244 Sys_Layout(ccomm, bs, &n, &N) 245 if na < (n+ng*b): raise ValueError( 246 "ghosts size %d, array size %d, and " 247 "vector local size %d block size %d" % 248 (toInt(ng), toInt(na), toInt(n), toInt(b))) 249 cdef PetscVec newvec = NULL 250 if bs == PETSC_DECIDE: 251 CHKERR( VecCreateGhostWithArray( 252 ccomm, n, N, ng, ig, sa, &newvec) ) 253 else: 254 CHKERR( VecCreateGhostBlockWithArray( 255 ccomm, bs, n, N, ng, ig, sa, &newvec) ) 256 PetscCLEAR(self.obj); self.vec = newvec 257 self.set_attr('__array__', array) 258 return self 259 260 def createShared(self, size, bsize=None, comm=None): 261 cdef MPI_Comm ccomm = def_Comm(comm, PETSC_COMM_DEFAULT) 262 cdef PetscInt bs=0, n=0, N=0 263 Vec_Sizes(size, bsize, &bs, &n, &N) 264 Sys_Layout(ccomm, bs, &n, &N) 265 cdef PetscVec newvec = NULL 266 CHKERR( VecCreateShared(ccomm, n, N, &newvec) ) 267 PetscCLEAR(self.obj); self.vec = newvec 268 if bs != PETSC_DECIDE: 269 CHKERR( VecSetBlockSize(self.vec, bs) ) 270 return self 271 272 def createNest(self, vecs, isets=None, comm=None): 273 vecs = list(vecs) 274 if isets: 275 isets = list(isets) 276 assert len(isets) == len(vecs) 277 else: 278 isets = None 279 cdef MPI_Comm ccomm = def_Comm(comm, PETSC_COMM_DEFAULT) 280 cdef Py_ssize_t i, m = len(vecs) 281 cdef PetscInt n = <PetscInt>m 282 cdef PetscVec *cvecs = NULL 283 cdef PetscIS *cisets = NULL 284 cdef object tmp1, tmp2 285 tmp1 = oarray_p(empty_p(n), NULL, <void**>&cvecs) 286 for i from 0 <= i < m: cvecs[i] = (<Vec?>vecs[i]).vec 287 if isets is not None: 288 tmp2 = oarray_p(empty_p(n), NULL, <void**>&cisets) 289 for i from 0 <= i < m: cisets[i] = (<IS?>isets[i]).iset 290 cdef PetscVec newvec = NULL 291 CHKERR( VecCreateNest(ccomm, n, cisets, cvecs,&newvec) ) 292 PetscCLEAR(self.obj); self.vec = newvec 293 return self 294 295 # 296 297 def setOptionsPrefix(self, prefix): 298 cdef const char *cval = NULL 299 prefix = str2bytes(prefix, &cval) 300 CHKERR( VecSetOptionsPrefix(self.vec, cval) ) 301 302 def getOptionsPrefix(self): 303 cdef const char *cval = NULL 304 CHKERR( VecGetOptionsPrefix(self.vec, &cval) ) 305 return bytes2str(cval) 306 307 def setFromOptions(self): 308 CHKERR( VecSetFromOptions(self.vec) ) 309 310 def setUp(self): 311 CHKERR( VecSetUp(self.vec) ) 312 return self 313 314 def setOption(self, option, flag): 315 CHKERR( VecSetOption(self.vec, option, flag) ) 316 317 def getType(self): 318 cdef PetscVecType cval = NULL 319 CHKERR( VecGetType(self.vec, &cval) ) 320 return bytes2str(cval) 321 322 def getSize(self): 323 cdef PetscInt N = 0 324 CHKERR( VecGetSize(self.vec, &N) ) 325 return toInt(N) 326 327 def getLocalSize(self): 328 cdef PetscInt n = 0 329 CHKERR( VecGetLocalSize(self.vec, &n) ) 330 return toInt(n) 331 332 def getSizes(self): 333 cdef PetscInt n = 0, N = 0 334 CHKERR( VecGetLocalSize(self.vec, &n) ) 335 CHKERR( VecGetSize(self.vec, &N) ) 336 return (toInt(n), toInt(N)) 337 338 def setBlockSize(self, bsize): 339 cdef PetscInt bs = asInt(bsize) 340 CHKERR( VecSetBlockSize(self.vec, bs) ) 341 342 def getBlockSize(self): 343 cdef PetscInt bs=0 344 CHKERR( VecGetBlockSize(self.vec, &bs) ) 345 return toInt(bs) 346 347 def getOwnershipRange(self): 348 cdef PetscInt low=0, high=0 349 CHKERR( VecGetOwnershipRange(self.vec, &low, &high) ) 350 return (toInt(low), toInt(high)) 351 352 def getOwnershipRanges(self): 353 cdef const PetscInt *rng = NULL 354 CHKERR( VecGetOwnershipRanges(self.vec, &rng) ) 355 cdef MPI_Comm comm = MPI_COMM_NULL 356 CHKERR( PetscObjectGetComm(<PetscObject>self.vec, &comm) ) 357 cdef int size = -1 358 CHKERR( MPI_Comm_size(comm, &size) ) 359 return array_i(size+1, rng) 360 361 def getBuffer(self, readonly=False): 362 if readonly: 363 return vec_getbuffer_r(self) 364 else: 365 return vec_getbuffer_w(self) 366 367 def getArray(self, readonly=False): 368 if readonly: 369 return vec_getarray_r(self) 370 else: 371 return vec_getarray_w(self) 372 373 def setArray(self, array): 374 vec_setarray(self, array) 375 376 def placeArray(self, array): 377 cdef PetscInt nv=0 378 cdef PetscInt na=0 379 cdef PetscScalar *a = NULL 380 CHKERR( VecGetLocalSize(self.vec, &nv) ) 381 array = oarray_s(array, &na, &a) 382 if (na != nv): raise ValueError( 383 "cannot place input array size %d, vector size %d" % 384 (toInt(na), toInt(nv))) 385 CHKERR( VecPlaceArray(self.vec, a) ) 386 self.set_attr('__placed_array__', array) 387 388 def resetArray(self, force=False): 389 cdef object array = None 390 array = self.get_attr('__placed_array__') 391 if array is None and not force: return None 392 CHKERR( VecResetArray(self.vec) ) 393 self.set_attr('__placed_array__', None) 394 return array 395 396 def getCUDAHandle(self, mode='rw'): 397 cdef PetscScalar *hdl = NULL 398 cdef const char *m = NULL 399 if mode is not None: mode = str2bytes(mode, &m) 400 if m == NULL or (m[0] == c'r' and m[1] == c'w'): 401 CHKERR( VecCUDAGetArray(self.vec, &hdl) ) 402 elif m[0] == c'r': 403 CHKERR( VecCUDAGetArrayRead(self.vec, <const PetscScalar**>&hdl) ) 404 elif m[0] == c'w': 405 CHKERR( VecCUDAGetArrayWrite(self.vec, &hdl) ) 406 else: 407 raise ValueError("Invalid mode: expected 'rw', 'r', or 'w'") 408 return <Py_uintptr_t>hdl 409 410 def restoreCUDAHandle(self, handle, mode='rw'): 411 cdef PetscScalar *hdl = <PetscScalar*>(<Py_uintptr_t>handle) 412 cdef const char *m = NULL 413 if mode is not None: mode = str2bytes(mode, &m) 414 if m == NULL or (m[0] == c'r' and m[1] == c'w'): 415 CHKERR( VecCUDARestoreArray(self.vec, &hdl) ) 416 elif m[0] == c'r': 417 CHKERR( VecCUDARestoreArrayRead(self.vec, <const PetscScalar**>&hdl) ) 418 elif m[0] == c'w': 419 CHKERR( VecCUDARestoreArrayWrite(self.vec, &hdl) ) 420 else: 421 raise ValueError("Invalid mode: expected 'rw', 'r', or 'w'") 422 423 def duplicate(self, array=None): 424 cdef Vec vec = type(self)() 425 CHKERR( VecDuplicate(self.vec, &vec.vec) ) 426 if array is not None: 427 vec_setarray(vec, array) 428 return vec 429 430 def copy(self, Vec result=None): 431 if result is None: 432 result = type(self)() 433 if result.vec == NULL: 434 CHKERR( VecDuplicate(self.vec, &result.vec) ) 435 CHKERR( VecCopy(self.vec, result.vec) ) 436 return result 437 438 def chop(self, tol): 439 cdef PetscReal rval = asReal(tol) 440 CHKERR( VecChop(self.vec, rval) ) 441 442 def load(self, Viewer viewer): 443 cdef MPI_Comm comm = MPI_COMM_NULL 444 cdef PetscObject obj = <PetscObject>(viewer.vwr) 445 if self.vec == NULL: 446 CHKERR( PetscObjectGetComm(obj, &comm) ) 447 CHKERR( VecCreate(comm, &self.vec) ) 448 CHKERR( VecLoad(self.vec, viewer.vwr) ) 449 return self 450 451 def equal(self, Vec vec): 452 cdef PetscBool flag = PETSC_FALSE 453 CHKERR( VecEqual(self.vec, vec.vec, &flag) ) 454 return toBool(flag) 455 456 def dot(self, Vec vec): 457 cdef PetscScalar sval = 0 458 CHKERR( VecDot(self.vec, vec.vec, &sval) ) 459 return toScalar(sval) 460 461 def dotBegin(self, Vec vec): 462 cdef PetscScalar sval = 0 463 CHKERR( VecDotBegin(self.vec, vec.vec, &sval) ) 464 465 def dotEnd(self, Vec vec): 466 cdef PetscScalar sval = 0 467 CHKERR( VecDotEnd(self.vec, vec.vec, &sval) ) 468 return toScalar(sval) 469 470 def tDot(self, Vec vec): 471 cdef PetscScalar sval = 0 472 CHKERR( VecTDot(self.vec, vec.vec, &sval) ) 473 return toScalar(sval) 474 475 def tDotBegin(self, Vec vec): 476 cdef PetscScalar sval = 0 477 CHKERR( VecTDotBegin(self.vec, vec.vec, &sval) ) 478 479 def tDotEnd(self, Vec vec): 480 cdef PetscScalar sval = 0 481 CHKERR( VecTDotEnd(self.vec, vec.vec, &sval) ) 482 return toScalar(sval) 483 484 def mDot(self, vecs, out=None): 485 <void>self; <void>vecs; <void>out; # unused 486 raise NotImplementedError 487 488 def mDotBegin(self, vecs, out=None): 489 <void>self; <void>vecs; <void>out; # unused 490 raise NotImplementedError 491 492 def mDotEnd(self, vecs, out=None): 493 <void>self; <void>vecs; <void>out; # unused 494 raise NotImplementedError 495 496 def mtDot(self, vecs, out=None): 497 <void>self; <void>vecs; <void>out; # unused 498 raise NotImplementedError 499 500 def mtDotBegin(self, vecs, out=None): 501 <void>self; <void>vecs; <void>out; # unused 502 raise NotImplementedError 503 504 def mtDotEnd(self, vecs, out=None): 505 <void>self; <void>vecs; <void>out; # unused 506 raise NotImplementedError 507 508 def norm(self, norm_type=None): 509 cdef PetscNormType norm_1_2 = PETSC_NORM_1_AND_2 510 cdef PetscNormType ntype = PETSC_NORM_2 511 if norm_type is not None: ntype = norm_type 512 cdef PetscReal rval[2] 513 CHKERR( VecNorm(self.vec, ntype, rval) ) 514 if ntype != norm_1_2: return toReal(rval[0]) 515 else: return (toReal(rval[0]), toReal(rval[1])) 516 517 def normBegin(self, norm_type=None): 518 cdef PetscNormType ntype = PETSC_NORM_2 519 if norm_type is not None: ntype = norm_type 520 cdef PetscReal dummy[2] 521 CHKERR( VecNormBegin(self.vec, ntype, dummy) ) 522 523 def normEnd(self, norm_type=None): 524 cdef PetscNormType norm_1_2 = PETSC_NORM_1_AND_2 525 cdef PetscNormType ntype = PETSC_NORM_2 526 if norm_type is not None: ntype = norm_type 527 cdef PetscReal rval[2] 528 CHKERR( VecNormEnd(self.vec, ntype, rval) ) 529 if ntype != norm_1_2: return toReal(rval[0]) 530 else: return (toReal(rval[0]), toReal(rval[1])) 531 532 def sum(self): 533 cdef PetscScalar sval = 0 534 CHKERR( VecSum(self.vec, &sval) ) 535 return toScalar(sval) 536 537 def min(self): 538 cdef PetscInt ival = 0 539 cdef PetscReal rval = 0 540 CHKERR( VecMin(self.vec, &ival, &rval) ) 541 return (toInt(ival), toReal(rval)) 542 543 def max(self): 544 cdef PetscInt ival = 0 545 cdef PetscReal rval = 0 546 CHKERR( VecMax(self.vec, &ival, &rval) ) 547 return (toInt(ival), toReal(rval)) 548 549 def normalize(self): 550 cdef PetscReal rval = 0 551 CHKERR( VecNormalize(self.vec, &rval) ) 552 return toReal(rval) 553 554 def reciprocal(self): 555 CHKERR( VecReciprocal(self.vec) ) 556 557 def exp(self): 558 CHKERR( VecExp(self.vec) ) 559 560 def log(self): 561 CHKERR( VecLog(self.vec) ) 562 563 def sqrtabs(self): 564 CHKERR( VecSqrtAbs(self.vec) ) 565 566 def abs(self): 567 CHKERR( VecAbs(self.vec) ) 568 569 def conjugate(self): 570 CHKERR( VecConjugate(self.vec) ) 571 572 def setRandom(self, Random random=None): 573 cdef PetscRandom rnd = NULL 574 if random is not None: rnd = random.rnd 575 CHKERR( VecSetRandom(self.vec, rnd) ) 576 577 def permute(self, IS order, invert=False): 578 cdef PetscBool cinvert = PETSC_FALSE 579 if invert: cinvert = PETSC_TRUE 580 CHKERR( VecPermute(self.vec, order.iset, cinvert) ) 581 582 def zeroEntries(self): 583 CHKERR( VecZeroEntries(self.vec) ) 584 585 def set(self, alpha): 586 cdef PetscScalar sval = asScalar(alpha) 587 CHKERR( VecSet(self.vec, sval) ) 588 589 def isset(self, IS idx, alpha): 590 cdef PetscScalar aval = asScalar(alpha) 591 CHKERR( VecISSet(self.vec, idx.iset, aval) ) 592 593 def scale(self, alpha): 594 cdef PetscScalar sval = asScalar(alpha) 595 CHKERR( VecScale(self.vec, sval) ) 596 597 def shift(self, alpha): 598 cdef PetscScalar sval = asScalar(alpha) 599 CHKERR( VecShift(self.vec, sval) ) 600 601 def chop(self, tol): 602 cdef PetscReal rval = asReal(tol) 603 CHKERR( VecChop(self.vec, rval) ) 604 605 def swap(self, Vec vec): 606 CHKERR( VecSwap(self.vec, vec.vec) ) 607 608 def axpy(self, alpha, Vec x): 609 cdef PetscScalar sval = asScalar(alpha) 610 CHKERR( VecAXPY(self.vec, sval, x.vec) ) 611 612 def isaxpy(self, IS idx, alpha, Vec x): 613 cdef PetscScalar sval = asScalar(alpha) 614 CHKERR( VecISAXPY(self.vec, idx.iset, sval, x.vec) ) 615 616 def aypx(self, alpha, Vec x): 617 cdef PetscScalar sval = asScalar(alpha) 618 CHKERR( VecAYPX(self.vec, sval, x.vec) ) 619 620 def axpby(self, alpha, beta, Vec y): 621 cdef PetscScalar sval1 = asScalar(alpha) 622 cdef PetscScalar sval2 = asScalar(beta) 623 CHKERR( VecAXPBY(self.vec, sval1, sval2, y.vec) ) 624 625 def waxpy(self, alpha, Vec x, Vec y): 626 cdef PetscScalar sval = asScalar(alpha) 627 CHKERR( VecWAXPY(self.vec, sval, x.vec, y.vec) ) 628 629 def maxpy(self, alphas, vecs): 630 cdef PetscInt n = 0 631 cdef PetscScalar *a = NULL 632 cdef PetscVec *v = NULL 633 cdef object tmp1 = iarray_s(alphas, &n, &a) 634 cdef object tmp2 = oarray_p(empty_p(n),NULL, <void**>&v) 635 assert n == len(vecs) 636 cdef Py_ssize_t i=0 637 for i from 0 <= i < n: 638 v[i] = (<Vec?>(vecs[i])).vec 639 CHKERR( VecMAXPY(self.vec, n, a, v) ) 640 641 def pointwiseMult(self, Vec x, Vec y): 642 CHKERR( VecPointwiseMult(self.vec, x.vec, y.vec) ) 643 644 def pointwiseDivide(self, Vec x, Vec y): 645 CHKERR( VecPointwiseDivide(self.vec, x.vec, y.vec) ) 646 647 def pointwiseMin(self, Vec x, Vec y): 648 CHKERR( VecPointwiseMin(self.vec, x.vec, y.vec) ) 649 650 def pointwiseMax(self, Vec x, Vec y): 651 CHKERR( VecPointwiseMax(self.vec, x.vec, y.vec) ) 652 653 def pointwiseMaxAbs(self, Vec x, Vec y): 654 CHKERR( VecPointwiseMaxAbs(self.vec, x.vec, y.vec) ) 655 656 def maxPointwiseDivide(self, Vec vec): 657 cdef PetscReal rval = 0 658 CHKERR( VecMaxPointwiseDivide(self.vec, vec.vec, &rval) ) 659 return toReal(rval) 660 661 def getValue(self, index): 662 cdef PetscInt ival = asInt(index) 663 cdef PetscScalar sval = 0 664 CHKERR( VecGetValues(self.vec, 1, &ival, &sval) ) 665 return toScalar(sval) 666 667 def getValues(self, indices, values=None): 668 return vecgetvalues(self.vec, indices, values) 669 670 def getValuesStagStencil(self, indices, values=None): 671 raise NotImplementedError('getValuesStagStencil not yet implemented in petsc4py') 672 673 def setValue(self, index, value, addv=None): 674 cdef PetscInt ival = asInt(index) 675 cdef PetscScalar sval = asScalar(value) 676 cdef PetscInsertMode caddv = insertmode(addv) 677 CHKERR( VecSetValues(self.vec, 1, &ival, &sval, caddv) ) 678 679 def setValues(self, indices, values, addv=None): 680 vecsetvalues(self.vec, indices, values, addv, 0, 0) 681 682 def setValuesBlocked(self, indices, values, addv=None): 683 vecsetvalues(self.vec, indices, values, addv, 1, 0) 684 685 def setValuesStagStencil(self, indices, values, addv=None): 686 raise NotImplementedError('setValuesStagStencil not yet implemented in petsc4py') 687 688 def setLGMap(self, LGMap lgmap): 689 CHKERR( VecSetLocalToGlobalMapping(self.vec, lgmap.lgm) ) 690 691 def getLGMap(self): 692 cdef LGMap cmap = LGMap() 693 CHKERR( VecGetLocalToGlobalMapping(self.vec, &cmap.lgm) ) 694 PetscINCREF(cmap.obj) 695 return cmap 696 697 def setValueLocal(self, index, value, addv=None): 698 cdef PetscInt ival = asInt(index) 699 cdef PetscScalar sval = asScalar(value) 700 cdef PetscInsertMode caddv = insertmode(addv) 701 CHKERR( VecSetValuesLocal(self.vec, 1, &ival, &sval, caddv) ) 702 703 def setValuesLocal(self, indices, values, addv=None): 704 vecsetvalues(self.vec, indices, values, addv, 0, 1) 705 706 def setValuesBlockedLocal(self, indices, values, addv=None): 707 vecsetvalues(self.vec, indices, values, addv, 1, 1) 708 709 def assemblyBegin(self): 710 CHKERR( VecAssemblyBegin(self.vec) ) 711 712 def assemblyEnd(self): 713 CHKERR( VecAssemblyEnd(self.vec) ) 714 715 def assemble(self): 716 CHKERR( VecAssemblyBegin(self.vec) ) 717 CHKERR( VecAssemblyEnd(self.vec) ) 718 719 # --- methods for strided vectors --- 720 721 def strideScale(self, field, alpha): 722 cdef PetscInt ival = asInt(field) 723 cdef PetscScalar sval = asScalar(alpha) 724 CHKERR( VecStrideScale(self.vec, ival, sval) ) 725 726 def strideSum(self, field): 727 cdef PetscInt ival = asInt(field) 728 cdef PetscScalar sval = 0 729 CHKERR( VecStrideSum(self.vec, ival, &sval) ) 730 return toScalar(sval) 731 732 def strideMin(self, field): 733 cdef PetscInt ival1 = asInt(field) 734 cdef PetscInt ival2 = 0 735 cdef PetscReal rval = 0 736 CHKERR( VecStrideMin(self.vec, ival1, &ival2, &rval) ) 737 return (toInt(ival2), toReal(rval)) 738 739 def strideMax(self, field): 740 cdef PetscInt ival1 = asInt(field) 741 cdef PetscInt ival2 = 0 742 cdef PetscReal rval = 0 743 CHKERR( VecStrideMax(self.vec, ival1, &ival2, &rval) ) 744 return (toInt(ival2), toReal(rval)) 745 746 def strideNorm(self, field, norm_type=None): 747 cdef PetscInt ival = asInt(field) 748 cdef PetscNormType norm_1_2 = PETSC_NORM_1_AND_2 749 cdef PetscNormType ntype = PETSC_NORM_2 750 if norm_type is not None: ntype = norm_type 751 cdef PetscReal rval[2] 752 CHKERR( VecStrideNorm(self.vec, ival, ntype, rval) ) 753 if ntype != norm_1_2: return toReal(rval[0]) 754 else: return (toReal(rval[0]), toReal(rval[1])) 755 756 def strideScatter(self, field, Vec vec, addv=None): 757 cdef PetscInt ival = asInt(field) 758 cdef PetscInsertMode caddv = insertmode(addv) 759 CHKERR( VecStrideScatter(self.vec, ival, vec.vec, caddv) ) 760 761 def strideGather(self, field, Vec vec, addv=None): 762 cdef PetscInt ival = asInt(field) 763 cdef PetscInsertMode caddv = insertmode(addv) 764 CHKERR( VecStrideGather(self.vec, ival, vec.vec, caddv) ) 765 766 # --- methods for vectors with ghost values --- 767 768 def localForm(self): 769 """ 770 Intended for use in context manager:: 771 772 with vec.localForm() as lf: 773 use(lf) 774 """ 775 return _Vec_LocalForm(self) 776 777 def ghostUpdateBegin(self, addv=None, mode=None): 778 cdef PetscInsertMode caddv = insertmode(addv) 779 cdef PetscScatterMode csctm = scattermode(mode) 780 CHKERR( VecGhostUpdateBegin(self.vec, caddv, csctm) ) 781 782 def ghostUpdateEnd(self, addv=None, mode=None): 783 cdef PetscInsertMode caddv = insertmode(addv) 784 cdef PetscScatterMode csctm = scattermode(mode) 785 CHKERR( VecGhostUpdateEnd(self.vec, caddv, csctm) ) 786 787 def ghostUpdate(self, addv=None, mode=None): 788 cdef PetscInsertMode caddv = insertmode(addv) 789 cdef PetscScatterMode csctm = scattermode(mode) 790 CHKERR( VecGhostUpdateBegin(self.vec, caddv, csctm) ) 791 CHKERR( VecGhostUpdateEnd(self.vec, caddv, csctm) ) 792 793 def setMPIGhost(self, ghosts): 794 "Alternative to createGhost()" 795 cdef PetscInt ng=0, *ig=NULL 796 ghosts = iarray_i(ghosts, &ng, &ig) 797 CHKERR( VecMPISetGhost(self.vec, ng, ig) ) 798 799 # 800 801 def getSubVector(self, IS iset, Vec subvec=None): 802 if subvec is None: subvec = Vec() 803 else: CHKERR( VecDestroy(&subvec.vec) ) 804 CHKERR( VecGetSubVector(self.vec, iset.iset, &subvec.vec) ) 805 return subvec 806 807 def restoreSubVector(self, IS iset, Vec subvec): 808 CHKERR( VecRestoreSubVector(self.vec, iset.iset, &subvec.vec) ) 809 810 def getNestSubVecs(self): 811 cdef PetscInt N=0 812 cdef PetscVec* sx=NULL 813 CHKERR( VecNestGetSubVecs(self.vec, &N, &sx) ) 814 output = [] 815 for i in range(N): 816 pyvec = Vec() 817 pyvec.vec = sx[i] 818 CHKERR( PetscObjectReference(<PetscObject> pyvec.vec) ) 819 output.append(pyvec) 820 821 return output 822 823 def setNestSubVecs(self, sx, idxm=None): 824 if idxm is None: idxm = range(len(sx)) 825 else: assert len(idxm) == len(sx) 826 cdef PetscInt N = 0 827 cdef PetscInt* cidxm = NULL 828 idxm = iarray_i(idxm, &N, &cidxm) 829 830 831 cdef PetscVec* csx = NULL 832 tmp = oarray_p(empty_p(N), NULL, <void**>&csx) 833 for i from 0 <= i < N: csx[i] = (<Vec?>sx[i]).vec 834 835 CHKERR( VecNestSetSubVecs(self.vec, N, cidxm, csx) ) 836 837 # 838 839 property sizes: 840 def __get__(self): 841 return self.getSizes() 842 def __set__(self, value): 843 self.setSizes(value) 844 845 property size: 846 def __get__(self): 847 return self.getSize() 848 849 property local_size: 850 def __get__(self): 851 return self.getLocalSize() 852 853 property block_size: 854 def __get__(self): 855 return self.getBlockSize() 856 857 property owner_range: 858 def __get__(self): 859 return self.getOwnershipRange() 860 861 property owner_ranges: 862 def __get__(self): 863 return self.getOwnershipRanges() 864 865 property buffer_w: 866 "Vec buffer (writable)" 867 def __get__(self): 868 return self.getBuffer() 869 870 property buffer_r: 871 "Vec buffer (read-only)" 872 def __get__(self): 873 return self.getBuffer(True) 874 875 property array_w: 876 "Vec array (writable)" 877 def __get__(self): 878 return self.getArray() 879 def __set__(self, value): 880 cdef buf = self.getBuffer() 881 with buf as array: array[:] = value 882 883 property array_r: 884 "Vec array (read-only)" 885 def __get__(self): 886 return self.getArray(True) 887 888 property buffer: 889 def __get__(self): 890 return self.buffer_w 891 892 property array: 893 def __get__(self): 894 return self.array_w 895 def __set__(self, value): 896 self.array_w = value 897 898 # --- NumPy array interface (legacy) --- 899 900 property __array_interface__: 901 def __get__(self): 902 cdef buf = self.getBuffer() 903 return buf.__array_interface__ 904 905# -------------------------------------------------------------------- 906 907del VecType 908del VecOption 909 910# -------------------------------------------------------------------- 911