Actual source code: rvector.c

  1: #define PETSCVEC_DLL

  3: /*
  4:      Provides the interface functions for vector operations that have PetscScalar/PetscReal in the signature
  5:    These are the vector functions the user calls.
  6: */
 7:  #include private/vecimpl.h

 11: /*@
 12:    VecMaxPointwiseDivide - Computes the maximum of the componentwise division max = max_i abs(x_i/y_i).

 14:    Collective on Vec

 16:    Input Parameters:
 17: .  x, y  - the vectors

 19:    Output Parameter:
 20: .  max - the result

 22:    Level: advanced

 24:    Notes: any subset of the x and may be the same vector.
 25:           if a particular y_i is zero, it is treated as 1 in the above formula

 27: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
 28: @*/
 29: PetscErrorCode  VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
 30: {

 40:   if (x->map.N != y->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
 41:   if (x->map.n != y->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

 43:   (*x->ops->maxpointwisedivide)(x,y,max);
 44:   return(0);
 45: }

 49: /*@
 50:    VecDot - Computes the vector dot product.

 52:    Collective on Vec

 54:    Input Parameters:
 55: .  x, y - the vectors

 57:    Output Parameter:
 58: .  val - the dot product

 60:    Performance Issues:
 61: +    per-processor memory bandwidth
 62: .    interprocessor latency
 63: -    work load inbalance that causes certain processes to arrive much earlier than
 64:      others

 66:    Notes for Users of Complex Numbers:
 67:    For complex vectors, VecDot() computes 
 68: $     val = (x,y) = y^H x,
 69:    where y^H denotes the conjugate transpose of y.

 71:    Use VecTDot() for the indefinite form
 72: $     val = (x,y) = y^T x,
 73:    where y^T denotes the transpose of y.

 75:    Level: intermediate

 77:    Concepts: inner product
 78:    Concepts: vector^inner product

 80: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd()
 81: @*/
 82: PetscErrorCode  VecDot(Vec x,Vec y,PetscScalar *val)
 83: {

 93:   if (x->map.N != y->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
 94:   if (x->map.n != y->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

 97:   (*x->ops->dot)(x,y,val);
 99:   return(0);
100: }

104: /*@
105:    VecNorm  - Computes the vector norm.

107:    Collective on Vec

109:    Input Parameters:
110: +  x - the vector
111: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
112:           NORM_1_AND_2, which computes both norms and stores them
113:           in a two element array.

115:    Output Parameter:
116: .  val - the norm 

118:    Notes:
119: $     NORM_1 denotes sum_i |x_i|
120: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
121: $     NORM_INFINITY denotes max_i |x_i|

123:    Level: intermediate

125:    Performance Issues:
126: +    per-processor memory bandwidth
127: .    interprocessor latency
128: -    work load inbalance that causes certain processes to arrive much earlier than
129:      others

131:    Compile Option:
132:    PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
133:  than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines 
134:  (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow. 

136:    Concepts: norm
137:    Concepts: vector^norm

139: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(), 
140:           VecNormBegin(), VecNormEnd()

142: @*/
143: PetscErrorCode  VecNorm(Vec x,NormType type,PetscReal *val)
144: {
145:   PetscTruth     flg;


153:   /*
154:    * Cached data?
155:    */
156:   if (type!=NORM_1_AND_2) {
157:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[type],*val,flg);
158:     if (flg) return(0);
159:   }

162:   (*x->ops->norm)(x,type,val);

165:   if (type!=NORM_1_AND_2) {
166:     PetscObjectComposedDataSetReal((PetscObject)x,NormIds[type],*val);
167:   }
168:   return(0);
169: }

173: /*@
174:    VecNormalize - Normalizes a vector by 2-norm. 

176:    Collective on Vec

178:    Input Parameters:
179: +  x - the vector

181:    Output Parameter:
182: .  x - the normalized vector
183: -  val - the vector norm before normalization

185:    Level: intermediate

187:    Concepts: vector^normalizing
188:    Concepts: normalizing^vector

190: @*/
191: PetscErrorCode  VecNormalize(Vec x,PetscReal *val)
192: {

200:   VecNorm(x,NORM_2,val);
201:   if (!*val) {
202:     PetscInfo(x,"Vector of zero norm can not be normalized; Returning only the zero norm\n");
203:   } else if (*val != 1.0) {
204:     PetscScalar tmp = 1.0/(*val);
205:     VecScale(x,tmp);
206:   }
208:   return(0);
209: }

213: /*@C
214:    VecMax - Determines the maximum vector component and its location.

216:    Collective on Vec

218:    Input Parameter:
219: .  x - the vector

221:    Output Parameters:
222: +  val - the maximum component
223: -  p - the location of val

225:    Notes:
226:    Returns the value PETSC_MIN and p = -1 if the vector is of length 0.

228:    Level: intermediate

230:    Concepts: maximum^of vector
231:    Concepts: vector^maximum value

233: .seealso: VecNorm(), VecMin()
234: @*/
235: PetscErrorCode  VecMax(Vec x,PetscInt *p,PetscReal *val)
236: {

244:   (*x->ops->max)(x,p,val);
246:   return(0);
247: }

251: /*@
252:    VecMin - Determines the minimum vector component and its location.

254:    Collective on Vec

256:    Input Parameters:
257: .  x - the vector

259:    Output Parameter:
260: +  val - the minimum component
261: -  p - the location of val

263:    Level: intermediate

265:    Notes:
266:    Returns the value PETSC_MAX and p = -1 if the vector is of length 0.

268:    Concepts: minimum^of vector
269:    Concepts: vector^minimum entry

271: .seealso: VecMax()
272: @*/
273: PetscErrorCode  VecMin(Vec x,PetscInt *p,PetscReal *val)
274: {

282:   (*x->ops->min)(x,p,val);
284:   return(0);
285: }

289: /*@
290:    VecTDot - Computes an indefinite vector dot product. That is, this
291:    routine does NOT use the complex conjugate.

293:    Collective on Vec

295:    Input Parameters:
296: .  x, y - the vectors

298:    Output Parameter:
299: .  val - the dot product

301:    Notes for Users of Complex Numbers:
302:    For complex vectors, VecTDot() computes the indefinite form
303: $     val = (x,y) = y^T x,
304:    where y^T denotes the transpose of y.

306:    Use VecDot() for the inner product
307: $     val = (x,y) = y^H x,
308:    where y^H denotes the conjugate transpose of y.

310:    Level: intermediate

312:    Concepts: inner product^non-Hermitian
313:    Concepts: vector^inner product
314:    Concepts: non-Hermitian inner product

316: .seealso: VecDot(), VecMTDot()
317: @*/
318: PetscErrorCode  VecTDot(Vec x,Vec y,PetscScalar *val)
319: {

329:   if (x->map.N != y->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
330:   if (x->map.n != y->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

333:   (*x->ops->tdot)(x,y,val);
335:   return(0);
336: }

340: /*@
341:    VecScale - Scales a vector. 

343:    Collective on Vec

345:    Input Parameters:
346: +  x - the vector
347: -  alpha - the scalar

349:    Output Parameter:
350: .  x - the scaled vector

352:    Note:
353:    For a vector with n components, VecScale() computes 
354: $      x[i] = alpha * x[i], for i=1,...,n.

356:    Level: intermediate

358:    Concepts: vector^scaling
359:    Concepts: scaling^vector

361: @*/
362: PetscErrorCode  VecScale (Vec x, PetscScalar alpha)
363: {
364:   PetscReal      norms[4] = {0.0,0.0,0.0, 0.0};
365:   PetscTruth     flgs[4];
367:   PetscInt       i;

372:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
374:   (*x->ops->scale)(x,alpha);

376:   /*
377:    * Update cached data
378:    */
379:   for (i=0; i<4; i++) {
380:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
381:   }

383:   /* in general we consider this object touched */
384:   PetscObjectStateIncrease((PetscObject)x);

386:   for (i=0; i<4; i++) {
387:     if (flgs[i]) {
388:       PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],PetscAbsScalar(alpha)*norms[i]);
389:     }
390:   }

393:   return(0);
394: }


399: /*@
400:    VecSet - Sets all components of a vector to a single scalar value. 

402:    Collective on Vec

404:    Input Parameters:
405: +  x  - the vector
406: -  alpha - the scalar

408:    Output Parameter:
409: .  x  - the vector

411:    Note:
412:    For a vector of dimension n, VecSet() computes
413: $     x[i] = alpha, for i=1,...,n,
414:    so that all vector entries then equal the identical
415:    scalar value, alpha.  Use the more general routine
416:    VecSetValues() to set different vector entries.

418:    You CANNOT call this after you have called VecSetValues() but before you call 
419:    VecAssemblyBegin/End().

421:    Level: beginner

423: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()

425:    Concepts: vector^setting to constant

427: @*/
428: PetscErrorCode  VecSet(Vec x,PetscScalar alpha)
429: {
430:   PetscReal      val;

436:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"You cannot call this after you have called VecSetValues() but\n before you have called VecAssemblyBegin/End()");
437: #if defined (PETSC_USE_DEBUG)
438:  {
439:    PetscReal alpha_local,alpha_max;
440:    alpha_local = PetscAbsScalar(alpha);
441:    MPI_Allreduce(&alpha_local,&alpha_max,1,MPIU_REAL,MPI_MAX,x->comm);
442:    if (alpha_local != alpha_max) SETERRQ(PETSC_ERR_ARG_WRONG,"Same value should be used across all processors");
443:  }
444: #endif

447:   (*x->ops->set)(x,alpha);

450:   /*
451:    * Update cached data
452:    */
453:   /* in general we consider this object touched */
454:   PetscObjectStateIncrease((PetscObject)x);

456:   /* however, norms can be simply set */
457:   val = PetscAbsScalar(alpha);
458:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_1],x->map.N * val);
459:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_INFINITY],val);
460:   val = sqrt((double)x->map.N) * val;
461:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_2],val);
462:   PetscObjectComposedDataSetReal((PetscObject)x,NormIds[NORM_FROBENIUS],val);
463:   return(0);
464: }


469: /*@
470:    VecAXPY - Computes y = alpha x + y. 

472:    Collective on Vec

474:    Input Parameters:
475: +  alpha - the scalar
476: -  x, y  - the vectors

478:    Output Parameter:
479: .  y - output vector

481:    Level: intermediate

483:    Concepts: vector^BLAS
484:    Concepts: BLAS

486: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
487: @*/
488: PetscErrorCode  VecAXPY(Vec y,PetscScalar alpha,Vec x)
489: {

498:   if (x->map.N != y->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
499:   if (x->map.n != y->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

502:   (*y->ops->axpy)(y,alpha,x);
504:   PetscObjectStateIncrease((PetscObject)y);
505:   return(0);
506: }

510: /*@
511:    VecAXPBY - Computes y = alpha x + beta y. 

513:    Collective on Vec

515:    Input Parameters:
516: +  alpha,beta - the scalars
517: -  x, y  - the vectors

519:    Output Parameter:
520: .  y - output vector

522:    Level: intermediate

524:    Concepts: BLAS
525:    Concepts: vector^BLAS

527: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
528: @*/
529: PetscErrorCode  VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
530: {

539:   if (x->map.N != y->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
540:   if (x->map.n != y->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

543:   (*y->ops->axpby)(y,alpha,beta,x);
545:   PetscObjectStateIncrease((PetscObject)y);
546:   return(0);
547: }

551: /*@
552:    VecAYPX - Computes y = x + alpha y.

554:    Collective on Vec

556:    Input Parameters:
557: +  alpha - the scalar
558: -  x, y  - the vectors

560:    Output Parameter:
561: .  y - output vector

563:    Level: intermediate

565:    Concepts: vector^BLAS
566:    Concepts: BLAS

568: .seealso: VecAXPY(), VecWAXPY()
569: @*/
570: PetscErrorCode  VecAYPX(Vec y,PetscScalar alpha,Vec x)
571: {

580:   if (x->map.N != y->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
581:   if (x->map.n != y->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

584:    (*y->ops->aypx)(y,alpha,x);
586:   PetscObjectStateIncrease((PetscObject)y);
587:   return(0);
588: }


593: /*@
594:    VecWAXPY - Computes w = alpha x + y.

596:    Collective on Vec

598:    Input Parameters:
599: +  alpha - the scalar
600: -  x, y  - the vectors

602:    Output Parameter:
603: .  w - the result

605:    Level: intermediate

607:    Concepts: vector^BLAS
608:    Concepts: BLAS

610: .seealso: VecAXPY(), VecAYPX(), VecAXPBY()
611: @*/
612: PetscErrorCode  VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
613: {

625:   if (x->map.N != y->map.N || x->map.N != w->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
626:   if (x->map.n != y->map.n || x->map.n != w->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

629:    (*w->ops->waxpy)(w,alpha,x,y);
631:   PetscObjectStateIncrease((PetscObject)w);
632:   return(0);
633: }


638: /*@
639:    VecSetValues - Inserts or adds values into certain locations of a vector. 

641:    Not Collective

643:    Input Parameters:
644: +  x - vector to insert in
645: .  ni - number of elements to add
646: .  ix - indices where to add
647: .  y - array of values
648: -  iora - either INSERT_VALUES or ADD_VALUES, where
649:    ADD_VALUES adds values to any existing entries, and
650:    INSERT_VALUES replaces existing entries with new values

652:    Notes: 
653:    VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.

655:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES 
656:    options cannot be mixed without intervening calls to the assembly
657:    routines.

659:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
660:    MUST be called after all calls to VecSetValues() have been completed.

662:    VecSetValues() uses 0-based indices in Fortran as well as in C.

664:    Negative indices may be passed in ix, these rows are 
665:    simply ignored. This allows easily inserting element load matrices
666:    with homogeneous Dirchlet boundary conditions that you don't want represented
667:    in the vector.

669:    Level: beginner

671:    Concepts: vector^setting values

673: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
674:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
675: @*/
676: PetscErrorCode  VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
677: {

686:   (*x->ops->setvalues)(x,ni,ix,y,iora);
688:   PetscObjectStateIncrease((PetscObject)x);
689:   return(0);
690: }

694: /*@
695:    VecGetValues - Gets values from certain locations of a vector. Currently 
696:           can only get values on the same processor

698:     Collective on Vec
699:  
700:    Input Parameters:
701: +  x - vector to get values from
702: .  ni - number of elements to get
703: -  ix - indices where to get them from (in global 1d numbering)

705:    Output Parameter:
706: .   y - array of values

708:    Notes: 
709:    The user provides the allocated array y; it is NOT allocated in this routine

711:    VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.

713:    VecAssemblyBegin() and VecAssemblyEnd()  MUST be called before calling this

715:    VecGetValues() uses 0-based indices in Fortran as well as in C.

717:    Level: beginner

719:    Concepts: vector^getting values

721: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecGetValuesLocal(),
722:            VecGetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecSetValues()
723: @*/
724: PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
725: {

733:   (*x->ops->getvalues)(x,ni,ix,y);
734:   return(0);
735: }

739: /*@
740:    VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector. 

742:    Not Collective

744:    Input Parameters:
745: +  x - vector to insert in
746: .  ni - number of blocks to add
747: .  ix - indices where to add in block count, rather than element count
748: .  y - array of values
749: -  iora - either INSERT_VALUES or ADD_VALUES, where
750:    ADD_VALUES adds values to any existing entries, and
751:    INSERT_VALUES replaces existing entries with new values

753:    Notes: 
754:    VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j], 
755:    for j=0,...,bs, for i=0,...,ni-1. where bs was set with VecSetBlockSize().

757:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES 
758:    options cannot be mixed without intervening calls to the assembly
759:    routines.

761:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
762:    MUST be called after all calls to VecSetValuesBlocked() have been completed.

764:    VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.

766:    Negative indices may be passed in ix, these rows are 
767:    simply ignored. This allows easily inserting element load matrices
768:    with homogeneous Dirchlet boundary conditions that you don't want represented
769:    in the vector.

771:    Level: intermediate

773:    Concepts: vector^setting values blocked

775: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
776:            VecSetValues()
777: @*/
778: PetscErrorCode  VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
779: {

788:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
790:   PetscObjectStateIncrease((PetscObject)x);
791:   return(0);
792: }


797: /*@
798:    VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
799:    using a local ordering of the nodes. 

801:    Not Collective

803:    Input Parameters:
804: +  x - vector to insert in
805: .  ni - number of elements to add
806: .  ix - indices where to add
807: .  y - array of values
808: -  iora - either INSERT_VALUES or ADD_VALUES, where
809:    ADD_VALUES adds values to any existing entries, and
810:    INSERT_VALUES replaces existing entries with new values

812:    Level: intermediate

814:    Notes: 
815:    VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.

817:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES 
818:    options cannot be mixed without intervening calls to the assembly
819:    routines.

821:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
822:    MUST be called after all calls to VecSetValuesLocal() have been completed.

824:    VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.

826:    Concepts: vector^setting values with local numbering

828: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
829:            VecSetValuesBlockedLocal()
830: @*/
831: PetscErrorCode  VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
832: {
834:   PetscInt       lixp[128],*lix = lixp;


843:   if (!x->ops->setvalueslocal) {
844:     if (!x->mapping) {
845:       SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
846:     }
847:     if (ni > 128) {
848:       PetscMalloc(ni*sizeof(PetscInt),&lix);
849:     }
850:     ISLocalToGlobalMappingApply(x->mapping,ni,(PetscInt*)ix,lix);
851:     (*x->ops->setvalues)(x,ni,lix,y,iora);
852:     if (ni > 128) {
853:       PetscFree(lix);
854:     }
855:   } else {
856:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
857:   }
859:   PetscObjectStateIncrease((PetscObject)x);
860:   return(0);
861: }

865: /*@
866:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
867:    using a local ordering of the nodes. 

869:    Not Collective

871:    Input Parameters:
872: +  x - vector to insert in
873: .  ni - number of blocks to add
874: .  ix - indices where to add in block count, not element count
875: .  y - array of values
876: -  iora - either INSERT_VALUES or ADD_VALUES, where
877:    ADD_VALUES adds values to any existing entries, and
878:    INSERT_VALUES replaces existing entries with new values

880:    Level: intermediate

882:    Notes: 
883:    VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j], 
884:    for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().

886:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES 
887:    options cannot be mixed without intervening calls to the assembly
888:    routines.

890:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd() 
891:    MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.

893:    VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.


896:    Concepts: vector^setting values blocked with local numbering

898: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(), 
899:            VecSetLocalToGlobalMappingBlock()
900: @*/
901: PetscErrorCode  VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
902: {
904:   PetscInt       lixp[128],*lix = lixp;

911:   if (!x->bmapping) {
912:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMappingBlock()");
913:   }
914:   if (ni > 128) {
915:     PetscMalloc(ni*sizeof(PetscInt),&lix);
916:   }

919:   ISLocalToGlobalMappingApply(x->bmapping,ni,(PetscInt*)ix,lix);
920:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
922:   if (ni > 128) {
923:     PetscFree(lix);
924:   }
925:   PetscObjectStateIncrease((PetscObject)x);
926:   return(0);
927: }



933: /*@
934:    VecMTDot - Computes indefinite vector multiple dot products. 
935:    That is, it does NOT use the complex conjugate.

937:    Collective on Vec

939:    Input Parameters:
940: +  x - one vector
941: .  nv - number of vectors
942: -  y - array of vectors.  Note that vectors are pointers

944:    Output Parameter:
945: .  val - array of the dot products

947:    Notes for Users of Complex Numbers:
948:    For complex vectors, VecMTDot() computes the indefinite form
949: $      val = (x,y) = y^T x,
950:    where y^T denotes the transpose of y.

952:    Use VecMDot() for the inner product
953: $      val = (x,y) = y^H x,
954:    where y^H denotes the conjugate transpose of y.

956:    Level: intermediate

958:    Concepts: inner product^multiple
959:    Concepts: vector^multiple inner products

961: .seealso: VecMDot(), VecTDot()
962: @*/
963: PetscErrorCode  VecMTDot(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
964: {

975:   if (x->map.N != (*y)->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
976:   if (x->map.n != (*y)->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

979:   (*x->ops->mtdot)(x,nv,y,val);
981:   return(0);
982: }

986: /*@
987:    VecMDot - Computes vector multiple dot products. 

989:    Collective on Vec

991:    Input Parameters:
992: +  x - one vector
993: .  nv - number of vectors
994: -  y - array of vectors. 

996:    Output Parameter:
997: .  val - array of the dot products

999:    Notes for Users of Complex Numbers:
1000:    For complex vectors, VecMDot() computes 
1001: $     val = (x,y) = y^H x,
1002:    where y^H denotes the conjugate transpose of y.

1004:    Use VecMTDot() for the indefinite form
1005: $     val = (x,y) = y^T x,
1006:    where y^T denotes the transpose of y.

1008:    Level: intermediate

1010:    Concepts: inner product^multiple
1011:    Concepts: vector^multiple inner products

1013: .seealso: VecMTDot(), VecDot()
1014: @*/
1015: PetscErrorCode  VecMDot(Vec x,PetscInt nv,const Vec y[],PetscScalar *val)
1016: {

1027:   if (x->map.N != (*y)->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1028:   if (x->map.n != (*y)->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1031:   (*x->ops->mdot)(x,nv,y,val);
1033:   return(0);
1034: }

1038: /*@
1039:    VecMAXPY - Computes y = y + sum alpha[j] x[j]

1041:    Collective on Vec

1043:    Input Parameters:
1044: +  nv - number of scalars and x-vectors
1045: .  alpha - array of scalars
1046: .  y - one vector
1047: -  x - array of vectors

1049:    Level: intermediate

1051:    Concepts: BLAS

1053: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
1054: @*/
1055: PetscErrorCode  VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
1056: {

1067:   if (y->map.N != (*x)->map.N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1068:   if (y->map.n != (*x)->map.n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1071:   (*y->ops->maxpy)(y,nv,alpha,x);
1073:   PetscObjectStateIncrease((PetscObject)y);
1074:   return(0);
1075: }

1077: /*MC
1078:    VecGetArray - Returns a pointer to a contiguous array that contains this 
1079:    processor's portion of the vector data. For the standard PETSc
1080:    vectors, VecGetArray() returns a pointer to the local data array and
1081:    does not use any copies. If the underlying vector data is not stored
1082:    in a contiquous array this routine will copy the data to a contiquous
1083:    array and return a pointer to that. You MUST call VecRestoreArray() 
1084:    when you no longer need access to the array.

1086:    Synopsis:
1087:    PetscErrorCode VecGetArray(Vec x,PetscScalar *a[])

1089:    Not Collective

1091:    Input Parameter:
1092: .  x - the vector

1094:    Output Parameter:
1095: .  a - location to put pointer to the array

1097:    Fortran Note:
1098:    This routine is used differently from Fortran 77
1099: $    Vec         x
1100: $    PetscScalar x_array(1)
1101: $    PetscOffset i_x
1102: $    PetscErrorCode ierr
1103: $       call VecGetArray(x,x_array,i_x,ierr)
1104: $
1105: $   Access first local entry in vector with
1106: $      value = x_array(i_x + 1)
1107: $
1108: $      ...... other code
1109: $       call VecRestoreArray(x,x_array,i_x,ierr)
1110:    For Fortran 90 see VecGetArrayF90()

1112:    See the Fortran chapter of the users manual and 
1113:    petsc/src/snes/examples/tutorials/ex5f.F for details.

1115:    Level: beginner

1117:    Concepts: vector^accessing local values

1119: .seealso: VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(), VecGetArray2d()
1120: M*/
1123: PetscErrorCode VecGetArray_Private(Vec x,PetscScalar *a[])
1124: {

1131:   (*x->ops->getarray)(x,a);
1132:   return(0);
1133: }


1138: /*@C
1139:    VecGetArrays - Returns a pointer to the arrays in a set of vectors
1140:    that were created by a call to VecDuplicateVecs().  You MUST call
1141:    VecRestoreArrays() when you no longer need access to the array.

1143:    Not Collective

1145:    Input Parameter:
1146: +  x - the vectors
1147: -  n - the number of vectors

1149:    Output Parameter:
1150: .  a - location to put pointer to the array

1152:    Fortran Note:
1153:    This routine is not supported in Fortran.

1155:    Level: intermediate

1157: .seealso: VecGetArray(), VecRestoreArrays()
1158: @*/
1159: PetscErrorCode  VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1160: {
1162:   PetscInt       i;
1163:   PetscScalar    **q;

1169:   if (n <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
1170:   PetscMalloc(n*sizeof(PetscScalar*),&q);
1171:   for (i=0; i<n; ++i) {
1172:     VecGetArray(x[i],&q[i]);
1173:   }
1174:   *a = q;
1175:   return(0);
1176: }

1180: /*@C
1181:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
1182:    has been called.

1184:    Not Collective

1186:    Input Parameters:
1187: +  x - the vector
1188: .  n - the number of vectors
1189: -  a - location of pointer to arrays obtained from VecGetArrays()

1191:    Notes:
1192:    For regular PETSc vectors this routine does not involve any copies. For
1193:    any special vectors that do not store local vector data in a contiguous
1194:    array, this routine will copy the data back into the underlying 
1195:    vector data structure from the arrays obtained with VecGetArrays().

1197:    Fortran Note:
1198:    This routine is not supported in Fortran.

1200:    Level: intermediate

1202: .seealso: VecGetArrays(), VecRestoreArray()
1203: @*/
1204: PetscErrorCode  VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
1205: {
1207:   PetscInt       i;
1208:   PetscScalar    **q = *a;


1215:   for(i=0;i<n;++i) {
1216:     VecRestoreArray(x[i],&q[i]);
1217:  }
1218:   PetscFree(q);
1219:   return(0);
1220: }

1222: /*MC
1223:    VecRestoreArray - Restores a vector after VecGetArray() has been called.

1225:    Not Collective

1227:    Synopsis:
1228:    PetscErrorCode VecRestoreArray(Vec x,PetscScalar *a[])

1230:    Input Parameters:
1231: +  x - the vector
1232: -  a - location of pointer to array obtained from VecGetArray()

1234:    Level: beginner

1236:    Notes:
1237:    For regular PETSc vectors this routine does not involve any copies. For
1238:    any special vectors that do not store local vector data in a contiguous
1239:    array, this routine will copy the data back into the underlying 
1240:    vector data structure from the array obtained with VecGetArray().

1242:    This routine actually zeros out the a pointer. This is to prevent accidental
1243:    us of the array after it has been restored. If you pass null for a it will 
1244:    not zero the array pointer a.

1246:    Fortran Note:
1247:    This routine is used differently from Fortran 77
1248: $    Vec         x
1249: $    PetscScalar x_array(1)
1250: $    PetscOffset i_x
1251: $    PetscErrorCode ierr
1252: $       call VecGetArray(x,x_array,i_x,ierr)
1253: $
1254: $   Access first local entry in vector with
1255: $      value = x_array(i_x + 1)
1256: $
1257: $      ...... other code
1258: $       call VecRestoreArray(x,x_array,i_x,ierr)

1260:    See the Fortran chapter of the users manual and 
1261:    petsc/src/snes/examples/tutorials/ex5f.F for details.
1262:    For Fortran 90 see VecRestoreArrayF90()

1264: .seealso: VecGetArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(), VecRestoreArray2d()
1265: M*/
1268: PetscErrorCode VecRestoreArray_Private(Vec x,PetscScalar *a[])
1269: {

1276: #if defined(PETSC_USE_DEBUG)
1277:   CHKMEMQ;
1278: #endif
1279:   if (x->ops->restorearray) {
1280:     (*x->ops->restorearray)(x,a);
1281:   }
1282:   PetscObjectStateIncrease((PetscObject)x);
1283:   return(0);
1284: }


1289: /*@
1290:    VecPlaceArray - Allows one to replace the array in a vector with an
1291:    array provided by the user. This is useful to avoid copying an array
1292:    into a vector.

1294:    Not Collective

1296:    Input Parameters:
1297: +  vec - the vector
1298: -  array - the array

1300:    Notes:
1301:    You can return to the original array with a call to VecResetArray()

1303:    Level: developer

1305: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()

1307: @*/
1308: PetscErrorCode  VecPlaceArray(Vec vec,const PetscScalar array[])
1309: {

1316:   if (vec->ops->placearray) {
1317:     (*vec->ops->placearray)(vec,array);
1318:   } else {
1319:     SETERRQ(PETSC_ERR_SUP,"Cannot place array in this type of vector");
1320:   }
1321:   PetscObjectStateIncrease((PetscObject)vec);
1322:   return(0);
1323: }


1328: /*@C
1329:    VecReplaceArray - Allows one to replace the array in a vector with an
1330:    array provided by the user. This is useful to avoid copying an array
1331:    into a vector.

1333:    Not Collective

1335:    Input Parameters:
1336: +  vec - the vector
1337: -  array - the array

1339:    Notes:
1340:    This permanently replaces the array and frees the memory associated
1341:    with the old array.

1343:    The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
1344:    freed by the user. It will be freed when the vector is destroy. 

1346:    Not supported from Fortran

1348:    Level: developer

1350: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()

1352: @*/
1353: PetscErrorCode  VecReplaceArray(Vec vec,const PetscScalar array[])
1354: {

1360:   if (vec->ops->replacearray) {
1361:     (*vec->ops->replacearray)(vec,array);
1362:  } else {
1363:     SETERRQ(PETSC_ERR_SUP,"Cannot replace array in this type of vector");
1364:   }
1365:   PetscObjectStateIncrease((PetscObject)vec);
1366:   return(0);
1367: }

1369: /*MC
1370:     VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
1371:     and makes them accessible via a Fortran90 pointer.

1373:     Synopsis:
1374:     VecDuplicateVecsF90(Vec x,int n,{Vec, pointer :: y(:)},integer ierr)

1376:     Collective on Vec

1378:     Input Parameters:
1379: +   x - a vector to mimic
1380: -   n - the number of vectors to obtain

1382:     Output Parameters:
1383: +   y - Fortran90 pointer to the array of vectors
1384: -   ierr - error code

1386:     Example of Usage: 
1387: .vb
1388:     Vec x
1389:     Vec, pointer :: y(:)
1390:     ....
1391:     call VecDuplicateVecsF90(x,2,y,ierr)
1392:     call VecSet(y(2),alpha,ierr)
1393:     call VecSet(y(2),alpha,ierr)
1394:     ....
1395:     call VecDestroyVecsF90(y,2,ierr)
1396: .ve

1398:     Notes:
1399:     Not yet supported for all F90 compilers

1401:     Use VecDestroyVecsF90() to free the space.

1403:     Level: beginner

1405: .seealso:  VecDestroyVecsF90(), VecDuplicateVecs()

1407: M*/

1409: /*MC
1410:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
1411:     VecGetArrayF90().

1413:     Synopsis:
1414:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1416:     Not collective

1418:     Input Parameters:
1419: +   x - vector
1420: -   xx_v - the Fortran90 pointer to the array

1422:     Output Parameter:
1423: .   ierr - error code

1425:     Example of Usage: 
1426: .vb
1427:     PetscScalar, pointer :: xx_v(:)
1428:     ....
1429:     call VecGetArrayF90(x,xx_v,ierr)
1430:     a = xx_v(3)
1431:     call VecRestoreArrayF90(x,xx_v,ierr)
1432: .ve
1433:    
1434:     Notes:
1435:     Not yet supported for all F90 compilers

1437:     Level: beginner

1439: .seealso:  VecGetArrayF90(), VecGetArray(), VecRestoreArray()

1441: M*/

1443: /*MC
1444:     VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().

1446:     Synopsis:
1447:     VecDestroyVecsF90({Vec, pointer :: x(:)},integer n,integer ierr)

1449:     Input Parameters:
1450: +   x - pointer to array of vector pointers
1451: -   n - the number of vectors previously obtained

1453:     Output Parameter:
1454: .   ierr - error code

1456:     Notes:
1457:     Not yet supported for all F90 compilers

1459:     Level: beginner

1461: .seealso:  VecDestroyVecs(), VecDuplicateVecsF90()

1463: M*/

1465: /*MC
1466:     VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
1467:     vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
1468:     this routine is implementation dependent. You MUST call VecRestoreArrayF90() 
1469:     when you no longer need access to the array.

1471:     Synopsis:
1472:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

1474:     Not Collective 

1476:     Input Parameter:
1477: .   x - vector

1479:     Output Parameters:
1480: +   xx_v - the Fortran90 pointer to the array
1481: -   ierr - error code

1483:     Example of Usage: 
1484: .vb
1485:     PetscScalar, pointer :: xx_v(:)
1486:     ....
1487:     call VecGetArrayF90(x,xx_v,ierr)
1488:     a = xx_v(3)
1489:     call VecRestoreArrayF90(x,xx_v,ierr)
1490: .ve

1492:     Notes:
1493:     Not yet supported for all F90 compilers

1495:     Level: beginner

1497: .seealso:  VecRestoreArrayF90(), VecGetArray(), VecRestoreArray()

1499: M*/


1504: /*@C
1505:    VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this 
1506:    processor's portion of the vector data.  You MUST call VecRestoreArray2d() 
1507:    when you no longer need access to the array.

1509:    Not Collective

1511:    Input Parameter:
1512: +  x - the vector
1513: .  m - first dimension of two dimensional array
1514: .  n - second dimension of two dimensional array
1515: .  mstart - first index you will use in first coordinate direction (often 0)
1516: -  nstart - first index in the second coordinate direction (often 0)

1518:    Output Parameter:
1519: .  a - location to put pointer to the array

1521:    Level: developer

1523:   Notes:
1524:    For a vector obtained from DACreateLocalVector() mstart and nstart are likely
1525:    obtained from the corner indices obtained from DAGetGhostCorners() while for
1526:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). In both cases
1527:    the arguments from DAGet[Ghost]Corners() are reversed in the call to VecGetArray2d().
1528:    
1529:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

1531:    Concepts: vector^accessing local values as 2d array

1533: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1534:           VecRestoreArray2d(), DAVecGetArray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1535:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1536: @*/
1537: PetscErrorCode  VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1538: {
1540:   PetscInt       i,N;
1541:   PetscScalar    *aa;

1547:   VecGetLocalSize(x,&N);
1548:   if (m*n != N) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
1549:   VecGetArray(x,&aa);

1551:   PetscMalloc(m*sizeof(PetscScalar*),a);
1552:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
1553:   *a -= mstart;
1554:   return(0);
1555: }

1559: /*@C
1560:    VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.

1562:    Not Collective

1564:    Input Parameters:
1565: +  x - the vector
1566: .  m - first dimension of two dimensional array
1567: .  n - second dimension of the two dimensional array
1568: .  mstart - first index you will use in first coordinate direction (often 0)
1569: .  nstart - first index in the second coordinate direction (often 0)
1570: -  a - location of pointer to array obtained from VecGetArray2d()

1572:    Level: developer

1574:    Notes:
1575:    For regular PETSc vectors this routine does not involve any copies. For
1576:    any special vectors that do not store local vector data in a contiguous
1577:    array, this routine will copy the data back into the underlying 
1578:    vector data structure from the array obtained with VecGetArray().

1580:    This routine actually zeros out the a pointer. 

1582: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1583:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
1584:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1585: @*/
1586: PetscErrorCode  VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
1587: {
1589:   void           *dummy;

1595:   dummy = (void*)(*a + mstart);
1596:   PetscFree(dummy);
1597:   VecRestoreArray(x,PETSC_NULL);
1598:   return(0);
1599: }

1603: /*@C
1604:    VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this 
1605:    processor's portion of the vector data.  You MUST call VecRestoreArray1d() 
1606:    when you no longer need access to the array.

1608:    Not Collective

1610:    Input Parameter:
1611: +  x - the vector
1612: .  m - first dimension of two dimensional array
1613: -  mstart - first index you will use in first coordinate direction (often 0)

1615:    Output Parameter:
1616: .  a - location to put pointer to the array

1618:    Level: developer

1620:   Notes:
1621:    For a vector obtained from DACreateLocalVector() mstart are likely
1622:    obtained from the corner indices obtained from DAGetGhostCorners() while for
1623:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). 
1624:    
1625:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

1627: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1628:           VecRestoreArray2d(), DAVecGetArray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1629:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1630: @*/
1631: PetscErrorCode  VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
1632: {
1634:   PetscInt       N;

1640:   VecGetLocalSize(x,&N);
1641:   if (m != N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
1642:   VecGetArray(x,a);
1643:   *a  -= mstart;
1644:   return(0);
1645: }

1649: /*@C
1650:    VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.

1652:    Not Collective

1654:    Input Parameters:
1655: +  x - the vector
1656: .  m - first dimension of two dimensional array
1657: .  mstart - first index you will use in first coordinate direction (often 0)
1658: -  a - location of pointer to array obtained from VecGetArray21()

1660:    Level: developer

1662:    Notes:
1663:    For regular PETSc vectors this routine does not involve any copies. For
1664:    any special vectors that do not store local vector data in a contiguous
1665:    array, this routine will copy the data back into the underlying 
1666:    vector data structure from the array obtained with VecGetArray1d().

1668:    This routine actually zeros out the a pointer. 

1670:    Concepts: vector^accessing local values as 1d array

1672: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1673:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
1674:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
1675: @*/
1676: PetscErrorCode  VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
1677: {

1683:   VecRestoreArray(x,PETSC_NULL);
1684:   return(0);
1685: }


1690: /*@C
1691:    VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this 
1692:    processor's portion of the vector data.  You MUST call VecRestoreArray3d() 
1693:    when you no longer need access to the array.

1695:    Not Collective

1697:    Input Parameter:
1698: +  x - the vector
1699: .  m - first dimension of three dimensional array
1700: .  n - second dimension of three dimensional array
1701: .  p - third dimension of three dimensional array
1702: .  mstart - first index you will use in first coordinate direction (often 0)
1703: .  nstart - first index in the second coordinate direction (often 0)
1704: -  pstart - first index in the third coordinate direction (often 0)

1706:    Output Parameter:
1707: .  a - location to put pointer to the array

1709:    Level: developer

1711:   Notes:
1712:    For a vector obtained from DACreateLocalVector() mstart, nstart, and pstart are likely
1713:    obtained from the corner indices obtained from DAGetGhostCorners() while for
1714:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). In both cases
1715:    the arguments from DAGet[Ghost]Corners() are reversed in the call to VecGetArray3d().
1716:    
1717:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

1719:    Concepts: vector^accessing local values as 3d array

1721: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1722:           VecRestoreArray2d(), DAVecGetarray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1723:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1724: @*/
1725: PetscErrorCode  VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
1726: {
1728:   PetscInt       i,N,j;
1729:   PetscScalar    *aa,**b;

1735:   VecGetLocalSize(x,&N);
1736:   if (m*n*p != N) SETERRQ4(PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 3d array dimensions %D by %D by %D",N,m,n,p);
1737:   VecGetArray(x,&aa);

1739:   PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
1740:   b    = (PetscScalar **)((*a) + m);
1741:   for (i=0; i<m; i++)   (*a)[i] = b + i*n - nstart;
1742:   for (i=0; i<m; i++) {
1743:     for (j=0; j<n; j++) {
1744:       b[i*n+j] = aa + i*n*p + j*p - pstart;
1745:     }
1746:   }
1747:   *a -= mstart;
1748:   return(0);
1749: }

1753: /*@C
1754:    VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.

1756:    Not Collective

1758:    Input Parameters:
1759: +  x - the vector
1760: .  m - first dimension of three dimensional array
1761: .  n - second dimension of the three dimensional array
1762: .  p - third dimension of the three dimensional array
1763: .  mstart - first index you will use in first coordinate direction (often 0)
1764: .  nstart - first index in the second coordinate direction (often 0)
1765: .  pstart - first index in the third coordinate direction (often 0)
1766: -  a - location of pointer to array obtained from VecGetArray3d()

1768:    Level: developer

1770:    Notes:
1771:    For regular PETSc vectors this routine does not involve any copies. For
1772:    any special vectors that do not store local vector data in a contiguous
1773:    array, this routine will copy the data back into the underlying 
1774:    vector data structure from the array obtained with VecGetArray().

1776:    This routine actually zeros out the a pointer. 

1778: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1779:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
1780:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
1781: @*/
1782: PetscErrorCode  VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
1783: {
1785:   void           *dummy;

1791:   dummy = (void*)(*a + mstart);
1792:   PetscFree(dummy);
1793:   VecRestoreArray(x,PETSC_NULL);
1794:   return(0);
1795: }

1799: /*@C
1800:    VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this 
1801:    processor's portion of the vector data.  You MUST call VecRestoreArray4d() 
1802:    when you no longer need access to the array.

1804:    Not Collective

1806:    Input Parameter:
1807: +  x - the vector
1808: .  m - first dimension of four dimensional array
1809: .  n - second dimension of four dimensional array
1810: .  p - third dimension of four dimensional array
1811: .  q - fourth dimension of four dimensional array
1812: .  mstart - first index you will use in first coordinate direction (often 0)
1813: .  nstart - first index in the second coordinate direction (often 0)
1814: .  pstart - first index in the third coordinate direction (often 0)
1815: -  qstart - first index in the fourth coordinate direction (often 0)

1817:    Output Parameter:
1818: .  a - location to put pointer to the array

1820:    Level: beginner

1822:   Notes:
1823:    For a vector obtained from DACreateLocalVector() mstart, nstart, and pstart are likely
1824:    obtained from the corner indices obtained from DAGetGhostCorners() while for
1825:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). In both cases
1826:    the arguments from DAGet[Ghost}Corners() are reversed in the call to VecGetArray3d().
1827:    
1828:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

1830:    Concepts: vector^accessing local values as 3d array

1832: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
1833:           VecRestoreArray2d(), DAVecGetarray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
1834:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
1835: @*/
1836: PetscErrorCode  VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
1837: {
1839:   PetscInt       i,N,j,k;
1840:   PetscScalar    *aa,***b,**c;

1846:   VecGetLocalSize(x,&N);
1847:   if (m*n*p*q != N) SETERRQ5(PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 4d array dimensions %D by %D by %D by %D",N,m,n,p,q);
1848:   VecGetArray(x,&aa);

1850:   PetscMalloc(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p*sizeof(PetscScalar*),a);
1851:   b    = (PetscScalar ***)((*a) + m);
1852:   c    = (PetscScalar **)(b + m*n);
1853:   for (i=0; i<m; i++)   (*a)[i] = b + i*n - nstart;
1854:   for (i=0; i<m; i++) {
1855:     for (j=0; j<n; j++) {
1856:       b[i*n+j] = c + i*n*p + j*p - pstart;
1857:     }
1858:   }
1859:   for (i=0; i<m; i++) {
1860:     for (j=0; j<n; j++) {
1861:       for (k=0; k<p; k++) {
1862:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
1863:       }
1864:     }
1865:   }
1866:   *a -= mstart;
1867:   return(0);
1868: }

1872: /*@C
1873:    VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.

1875:    Not Collective

1877:    Input Parameters:
1878: +  x - the vector
1879: .  m - first dimension of four dimensional array
1880: .  n - second dimension of the four dimensional array
1881: .  p - third dimension of the four dimensional array
1882: .  q - fourth dimension of the four dimensional array
1883: .  mstart - first index you will use in first coordinate direction (often 0)
1884: .  nstart - first index in the second coordinate direction (often 0)
1885: .  pstart - first index in the third coordinate direction (often 0)
1886: .  qstart - first index in the fourth coordinate direction (often 0)
1887: -  a - location of pointer to array obtained from VecGetArray4d()

1889:    Level: beginner

1891:    Notes:
1892:    For regular PETSc vectors this routine does not involve any copies. For
1893:    any special vectors that do not store local vector data in a contiguous
1894:    array, this routine will copy the data back into the underlying 
1895:    vector data structure from the array obtained with VecGetArray().

1897:    This routine actually zeros out the a pointer. 

1899: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
1900:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
1901:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
1902: @*/
1903: PetscErrorCode  VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
1904: {
1906:   void           *dummy;

1912:   dummy = (void*)(*a + mstart);
1913:   PetscFree(dummy);
1914:   VecRestoreArray(x,PETSC_NULL);
1915:   return(0);
1916: }