Actual source code: bvec2.c

petsc-3.5.1 2014-07-24
Report Typos and Errors
  2: /*
  3:    Implements the sequential vectors.
  4: */

  6: #include <../src/vec/vec/impls/dvecimpl.h>          /*I "petscvec.h" I*/
  7: #include <../src/vec/vec/impls/mpi/pvecimpl.h> /* For VecView_MPI_HDF5 */
  8: #include <petscblaslapack.h>
  9: #include <petscthreadcomm.h>

 11: #if defined(PETSC_HAVE_HDF5)
 12: extern PetscErrorCode VecView_MPI_HDF5(Vec,PetscViewer);
 13: #endif

 15: #if defined(PETSC_THREADCOMM_ACTIVE)
 16: PetscErrorCode VecPointwiseMax_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
 17: {
 18:   PetscErrorCode    ierr;
 19:   PetscInt          *trstarts=win->map->trstarts;
 20:   PetscInt          i;
 21:   const PetscScalar *xx,*yy;
 22:   PetscScalar       *ww;

 24:   VecGetArrayRead(xin,&xx);
 25:   VecGetArrayRead(yin,&yy);
 26:   VecGetArray(win,&ww);

 28:   for (i=trstarts[thread_id]; i < trstarts[thread_id+1];i++) ww[i] = PetscMax(PetscRealPart(xx[i]),PetscRealPart(yy[i]));

 30:   VecRestoreArrayRead(xin,&xx);
 31:   VecRestoreArrayRead(yin,&yy);
 32:   VecRestoreArray(win,&ww);
 33:   return 0;
 34: }

 38: PetscErrorCode VecPointwiseMax_Seq(Vec win,Vec xin,Vec yin)
 39: {

 43:   PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseMax_kernel,win,xin,yin);
 44:   PetscLogFlops(win->map->n);
 45:   return(0);
 46: }
 47: #else
 50: PetscErrorCode VecPointwiseMax_Seq(Vec win,Vec xin,Vec yin)
 51: {
 53:   PetscInt       n = win->map->n,i;
 54:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

 57:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
 58:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
 59:   VecGetArray(win,&ww);

 61:   for (i=0; i<n; i++) ww[i] = PetscMax(PetscRealPart(xx[i]),PetscRealPart(yy[i]));

 63:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
 64:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
 65:   VecRestoreArray(win,&ww);
 66:   PetscLogFlops(n);
 67:   return(0);
 68: }
 69: #endif

 71: #if defined(PETSC_THREADCOMM_ACTIVE)
 72: PetscErrorCode VecPointwiseMin_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
 73: {
 74:   PetscErrorCode    ierr;
 75:   PetscInt          *trstarts=win->map->trstarts;
 76:   PetscInt          i;
 77:   const PetscScalar *xx,*yy;
 78:   PetscScalar       *ww;

 80:   VecGetArrayRead(xin,&xx);
 81:   VecGetArrayRead(yin,&yy);
 82:   VecGetArray(win,&ww);

 84:   for (i=trstarts[thread_id]; i < trstarts[thread_id+1];i++) ww[i] = PetscMin(PetscRealPart(xx[i]),PetscRealPart(yy[i]));

 86:   VecRestoreArrayRead(xin,&xx);
 87:   VecRestoreArrayRead(yin,&yy);
 88:   VecRestoreArray(win,&ww);
 89:   return 0;
 90: }

 94: PetscErrorCode VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)
 95: {

 99:   PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseMin_kernel,win,xin,yin);
100:   PetscLogFlops(win->map->n);
101:   return(0);
102: }
103: #else
106: PetscErrorCode VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)
107: {
109:   PetscInt       n = win->map->n,i;
110:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

113:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
114:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
115:   VecGetArray(win,&ww);

117:   for (i=0; i<n; i++) ww[i] = PetscMin(PetscRealPart(xx[i]),PetscRealPart(yy[i]));

119:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
120:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
121:   VecRestoreArray(win,&ww);
122:   PetscLogFlops(n);
123:   return(0);
124: }
125: #endif

127: #if defined(PETSC_THREADCOMM_ACTIVE)
128: PetscErrorCode VecPointwiseMaxAbs_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
129: {
130:   PetscErrorCode    ierr;
131:   PetscInt          *trstarts=win->map->trstarts;
132:   PetscInt          i;
133:   const PetscScalar *xx,*yy;
134:   PetscScalar       *ww;

136:   VecGetArrayRead(xin,&xx);
137:   VecGetArrayRead(yin,&yy);
138:   VecGetArray(win,&ww);

140:   for (i=trstarts[thread_id]; i < trstarts[thread_id+1];i++) ww[i] = PetscMax(PetscAbsScalar(xx[i]),PetscAbsScalar(yy[i]));

142:   VecRestoreArrayRead(xin,&xx);
143:   VecRestoreArrayRead(yin,&yy);
144:   VecRestoreArray(win,&ww);
145:   return 0;
146: }

150: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
151: {

155:   PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseMaxAbs_kernel,win,xin,yin);
156:   PetscLogFlops(win->map->n);
157:   return(0);
158: }
159: #else
162: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
163: {
165:   PetscInt       n = win->map->n,i;
166:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

169:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
170:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
171:   VecGetArray(win,&ww);

173:   for (i=0; i<n; i++) ww[i] = PetscMax(PetscAbsScalar(xx[i]),PetscAbsScalar(yy[i]));

175:   PetscLogFlops(n);
176:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
177:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
178:   VecRestoreArray(win,&ww);
179:   return(0);
180: }
181: #endif

183: #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>
184: #if defined(PETSC_THREADCOMM_ACTIVE)
185: PetscErrorCode VecPointwiseMult_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
186: {
187:   PetscErrorCode    ierr;
188:   PetscInt          *trstarts=win->map->trstarts;
189:   PetscInt          i;
190:   const PetscScalar *xx,*yy;
191:   PetscScalar       *ww;

193:   VecGetArrayRead(xin,&xx);
194:   VecGetArrayRead(yin,&yy);
195:   VecGetArray(win,&ww);
196:   if (ww == xx) {
197:     for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) ww[i] *= yy[i];
198:   } else if (ww == yy) {
199:     for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) ww[i] *= xx[i];
200:   } else {
201: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
202:     PetscInt start,n;
203:     start = trstarts[thread_id];
204:     n     = trstarts[thread_id+1] - trstarts[thread_id];
205:     fortranxtimesy_(xx+start,yy+start,ww+start,&n);
206:   }
207: #else
208:     for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) ww[i] = xx[i] * yy[i];
209:   }
210: #endif
211:   VecRestoreArrayRead(xin,&xx);
212:   VecRestoreArrayRead(yin,&yy);
213:   VecRestoreArray(win,&ww);
214:   return 0;
215: }

219: PetscErrorCode VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)
220: {

224:   PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseMult_kernel,win,xin,yin);
225:   PetscLogFlops(win->map->n);
226:   return(0);
227: }
228: #else
231: PetscErrorCode VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)
232: {
234:   PetscInt       n = win->map->n,i;
235:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

238:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
239:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
240:   VecGetArray(win,&ww);
241:   if (ww == xx) {
242:     for (i=0; i<n; i++) ww[i] *= yy[i];
243:   } else if (ww == yy) {
244:     for (i=0; i<n; i++) ww[i] *= xx[i];
245:   } else {
246: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
247:     fortranxtimesy_(xx,yy,ww,&n);
248: #else
249:     for (i=0; i<n; i++) ww[i] = xx[i] * yy[i];
250: #endif
251:   }
252:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
253:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
254:   VecRestoreArray(win,&ww);
255:   PetscLogFlops(n);
256:   return(0);
257: }
258: #endif

260: #if defined(PETSC_THREADCOMM_ACTIVE)
261: PetscErrorCode VecPointwiseDivide_kernel(PetscInt thread_id,Vec win,Vec xin,Vec yin)
262: {
263:   PetscErrorCode    ierr;
264:   PetscInt          *trstarts=win->map->trstarts;
265:   PetscInt          i;
266:   const PetscScalar *xx,*yy;
267:   PetscScalar       *ww;

269:   VecGetArrayRead(xin,&xx);
270:   VecGetArrayRead(yin,&yy);
271:   VecGetArray(win,&ww);

273:   for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) ww[i] = xx[i] / yy[i];

275:   VecRestoreArrayRead(xin,&xx);
276:   VecRestoreArrayRead(yin,&yy);
277:   VecRestoreArray(win,&ww);
278:   return 0;
279: }

283: PetscErrorCode VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)
284: {

288:   PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)win),(PetscThreadKernel)VecPointwiseDivide_kernel,win,xin,yin);
289:   PetscLogFlops(win->map->n);
290:   return(0);
291: }
292: #else
295: PetscErrorCode VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)
296: {
298:   PetscInt       n = win->map->n,i;
299:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

302:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
303:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
304:   VecGetArray(win,&ww);

306:   for (i=0; i<n; i++) ww[i] = xx[i] / yy[i];

308:   PetscLogFlops(n);
309:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
310:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
311:   VecRestoreArray(win,&ww);
312:   return(0);
313: }
314: #endif

316: #if defined(PETSC_THREADCOMM_ACTIVE)
317: PetscErrorCode VecSetRandom_kernel(PetscInt thread_id,Vec xin, PetscRandom r)
318: {
320:   PetscInt       *trstarts=xin->map->trstarts;
321:   PetscInt       i;
322:   PetscScalar    *xx;

324:   VecGetArray(xin,&xx);
325:   for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) {
326:     PetscRandomGetValue(r,&xx[i]);
327:   }
328:   VecRestoreArray(xin,&xx);
329:   return 0;
330: }

334: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
335: {

339:   PetscThreadCommRunKernel2(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecSetRandom_kernel,xin,r);
340:   return(0);
341: }
342: #else
345: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
346: {
348:   PetscInt       n = xin->map->n,i;
349:   PetscScalar    *xx;

352:   VecGetArray(xin,&xx);
353:   for (i=0; i<n; i++) {PetscRandomGetValue(r,&xx[i]);}
354:   VecRestoreArray(xin,&xx);
355:   return(0);
356: }
357: #endif

361: PetscErrorCode VecGetSize_Seq(Vec vin,PetscInt *size)
362: {
364:   *size = vin->map->n;
365:   return(0);
366: }


369: #if defined(PETSC_THREADCOMM_ACTIVE)
370: PetscErrorCode VecConjugate_kernel(PetscInt thread_id,Vec xin)
371: {
373:   PetscScalar    *x;
374:   PetscInt       *trstarts=xin->map->trstarts,i;

376:   VecGetArray(xin,&x);
377:   for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) x[i] = PetscConj(x[i]);
378:   VecRestoreArray(xin,&x);
379:   return 0;
380: }

384: PetscErrorCode VecConjugate_Seq(Vec xin)
385: {

389:   PetscThreadCommRunKernel1(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecConjugate_kernel,xin);
390:   return(0);
391: }
392: #else
395: PetscErrorCode VecConjugate_Seq(Vec xin)
396: {
397:   PetscScalar    *x;
398:   PetscInt       n = xin->map->n;

402:   VecGetArray(xin,&x);
403:   while (n-->0) {
404:     *x = PetscConj(*x);
405:     x++;
406:   }
407:   VecRestoreArray(xin,&x);
408:   return(0);
409: }
410: #endif

414: PetscErrorCode VecResetArray_Seq(Vec vin)
415: {
416:   Vec_Seq *v = (Vec_Seq*)vin->data;

419:   v->array         = v->unplacedarray;
420:   v->unplacedarray = 0;
421:   return(0);
422: }

424: #if defined(PETSC_THREADCOMM_ACTIVE)
425: PetscErrorCode VecCopy_kernel(PetscInt thread_id,Vec xin,Vec yin)
426: {
427:   PetscErrorCode    ierr;
428:   PetscInt          *trstarts=yin->map->trstarts;
429:   PetscScalar       *ya;
430:   const PetscScalar *xa;
431:   PetscInt          start,end,n;

433:   VecGetArrayRead(xin,&xa);
434:   VecGetArray(yin,&ya);
435:   start = trstarts[thread_id];
436:   end   = trstarts[thread_id+1];
437:   n     = end-start;
438:   PetscMemcpy(ya+start,xa+start,n*sizeof(PetscScalar));
439:   VecRestoreArrayRead(xin,&xa);
440:   VecRestoreArray(yin,&ya);
441:   return 0;
442: }

446: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
447: {

451:   if (xin != yin) {
452:     PetscThreadCommRunKernel2(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecCopy_kernel,xin,yin);
453:   }
454:   return(0);
455: }
456: #else
459: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
460: {
461:   PetscScalar       *ya;
462:   const PetscScalar *xa;
463:   PetscErrorCode    ierr;

466:   if (xin != yin) {
467:     VecGetArrayRead(xin,&xa);
468:     VecGetArray(yin,&ya);
469:     PetscMemcpy(ya,xa,xin->map->n*sizeof(PetscScalar));
470:     VecRestoreArrayRead(xin,&xa);
471:     VecRestoreArray(yin,&ya);
472:   }
473:   return(0);
474: }
475: #endif

477: #if defined(PETSC_THREADCOMM_ACTIVE)
478: PetscErrorCode VecSwap_kernel(PetscInt thread_id,Vec xin,Vec yin)
479: {
480:   PetscErrorCode   ierr;
481:   PetscInt         *trstarts=xin->map->trstarts;
482:   PetscInt         start,end,n;
483:   PetscBLASInt     one=1,bn;
484:   PetscScalar      *xa,*ya;

486:   VecGetArray(xin,&xa);
487:   VecGetArray(yin,&ya);
488:   start = trstarts[thread_id];
489:   end   = trstarts[thread_id+1];
490:   n     = end-start;
491:   PetscBLASIntCast(n,&bn);
492:   PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa+start,&one,ya+start,&one));
493:   VecRestoreArray(xin,&xa);
494:   VecRestoreArray(yin,&ya);
495:   return 0;
496: }

500: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
501: {

505:   if (xin != yin) {
506:     PetscThreadCommRunKernel2(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecSwap_kernel,xin,yin);
507:   }
508:   return(0);
509: }
510: #else

514: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
515: {
516:   PetscScalar    *ya, *xa;
518:   PetscBLASInt   one = 1,bn;

521:   if (xin != yin) {
522:     PetscBLASIntCast(xin->map->n,&bn);
523:     VecGetArray(xin,&xa);
524:     VecGetArray(yin,&ya);
525:     PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa,&one,ya,&one));
526:     VecRestoreArray(xin,&xa);
527:     VecRestoreArray(yin,&ya);
528:   }
529:   return(0);
530: }
531: #endif

533: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>
534: #if defined(PETSC_THREADCOMM_ACTIVE)
535: PetscErrorCode VecNorm_kernel(PetscInt thread_id,Vec xin,NormType* type_p,PetscThreadCommReduction red)
536: {
538:   PetscInt       *trstarts=xin->map->trstarts;
539:   PetscInt       start,end,n;
540:   PetscBLASInt   one = 1,bn;
541:   const PetscScalar *xx;
542:   PetscReal      z_loc;
543:   NormType       type=*type_p;

545:   start = trstarts[thread_id];
546:   end   = trstarts[thread_id+1];
547:   n     = end - start;
548:   PetscBLASIntCast(n,&bn);
549:   VecGetArrayRead(xin,&xx);
550:   if (type == NORM_2 || type == NORM_FROBENIUS) {
551:     z_loc = PetscRealPart(BLASdot_(&bn,xx+start,&one,xx+start,&one));
552:     PetscThreadReductionKernelPost(thread_id,red,(void*)&z_loc);
553:   } else if (type == NORM_INFINITY) {
554:     PetscInt i;
555:     PetscReal max=0.0,tmp;
556:     for (i=trstarts[thread_id]; i < trstarts[thread_id+1]; i++) {
557:       if ((tmp = PetscAbsScalar(xx[i])) > max) max = tmp;
558:       /* check special case of tmp == NaN */
559:       if (tmp != tmp) {max = tmp; break;}
560:     }
561:     PetscThreadReductionKernelPost(thread_id,red,(void*)&max);
562:   } else if (type == NORM_1) {
563:     PetscStackCallBLAS("BLASasum",z_loc = BLASasum_(&bn,xx+start,&one));
564:     PetscThreadReductionKernelPost(thread_id,red,(void*)&z_loc);
565:   }
566:   VecRestoreArrayRead(xin,&xx);
567:   return 0;
568: }

572: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
573: {
574:   PetscErrorCode           ierr;
575:   PetscThreadCommReduction red;

578:   if (type == NORM_2 || type == NORM_FROBENIUS) {
579:     PetscThreadReductionBegin(PetscObjectComm((PetscObject)xin),THREADCOMM_SUM,PETSC_REAL,1,&red);
580:     PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecNorm_kernel,xin,(void*)&type,red);
581:     PetscThreadReductionEnd(red,z);
582:     *z   = PetscSqrtReal(*z);
583:     PetscLogFlops(PetscMax(2*xin->map->n-1.0,0.0));
584:   } else if (type == NORM_INFINITY) {
585:     PetscThreadReductionBegin(PetscObjectComm((PetscObject)xin),THREADCOMM_MAX,PETSC_REAL,1,&red);
586:     PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecNorm_kernel,xin,(void*)&type,red);
587:     PetscThreadReductionEnd(red,z);
588:   } else if (type == NORM_1) {
589:     PetscThreadReductionBegin(PetscObjectComm((PetscObject)xin),THREADCOMM_SUM,PETSC_REAL,1,&red);
590:     PetscThreadCommRunKernel3(PetscObjectComm((PetscObject)xin),(PetscThreadKernel)VecNorm_kernel,xin,(void*)&type,red);
591:     PetscThreadReductionEnd(red,z);
592:     PetscLogFlops(PetscMax(xin->map->n-1.0,0.0));
593:   } else if (type == NORM_1_AND_2) {
594:     VecNorm_Seq(xin,NORM_1,z);
595:     VecNorm_Seq(xin,NORM_2,z+1);
596:   }
597:   return(0);
598: }
599: #else

603: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
604: {
605:   const PetscScalar *xx;
606:   PetscErrorCode    ierr;
607:   PetscInt          n = xin->map->n;
608:   PetscBLASInt      one = 1, bn;

611:   PetscBLASIntCast(n,&bn);
612:   if (type == NORM_2 || type == NORM_FROBENIUS) {
613:     VecGetArrayRead(xin,&xx);
614:     *z   = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
615:     *z   = PetscSqrtReal(*z);
616:     VecRestoreArrayRead(xin,&xx);
617:     PetscLogFlops(PetscMax(2.0*n-1,0.0));
618:   } else if (type == NORM_INFINITY) {
619:     PetscInt  i;
620:     PetscReal max = 0.0,tmp;

622:     VecGetArrayRead(xin,&xx);
623:     for (i=0; i<n; i++) {
624:       if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
625:       /* check special case of tmp == NaN */
626:       if (tmp != tmp) {max = tmp; break;}
627:       xx++;
628:     }
629:     VecRestoreArrayRead(xin,&xx);
630:     *z   = max;
631:   } else if (type == NORM_1) {
632:     VecGetArrayRead(xin,&xx);
633:     PetscStackCallBLAS("BLASasum",*z   = BLASasum_(&bn,xx,&one));
634:     VecRestoreArrayRead(xin,&xx);
635:     PetscLogFlops(PetscMax(n-1.0,0.0));
636:   } else if (type == NORM_1_AND_2) {
637:     VecNorm_Seq(xin,NORM_1,z);
638:     VecNorm_Seq(xin,NORM_2,z+1);
639:   }
640:   return(0);
641: }
642: #endif

646: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
647: {
648:   PetscErrorCode    ierr;
649:   PetscInt          i,n = xin->map->n;
650:   const char        *name;
651:   PetscViewerFormat format;
652:   const PetscScalar *xv;

655:   VecGetArrayRead(xin,&xv);
656:   PetscViewerGetFormat(viewer,&format);
657:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
658:     PetscObjectGetName((PetscObject)xin,&name);
659:     PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
660:     for (i=0; i<n; i++) {
661: #if defined(PETSC_USE_COMPLEX)
662:       if (PetscImaginaryPart(xv[i]) > 0.0) {
663:         PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
664:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
665:         PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
666:       } else {
667:         PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)PetscRealPart(xv[i]));
668:       }
669: #else
670:       PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double) xv[i]);
671: #endif
672:     }
673:     PetscViewerASCIIPrintf(viewer,"];\n");
674:   } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
675:     for (i=0; i<n; i++) {
676: #if defined(PETSC_USE_COMPLEX)
677:       PetscViewerASCIIPrintf(viewer,"%18.16e %18.16e\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
678: #else
679:       PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)xv[i]);
680: #endif
681:     }
682:   } else if (format == PETSC_VIEWER_ASCII_VTK || format == PETSC_VIEWER_ASCII_VTK_CELL) {
683:     /*
684:        state 0: No header has been output
685:        state 1: Only POINT_DATA has been output
686:        state 2: Only CELL_DATA has been output
687:        state 3: Output both, POINT_DATA last
688:        state 4: Output both, CELL_DATA last
689:     */
690:     static PetscInt stateId = -1;
691:     int outputState = 0;
692:     PetscBool  hasState;
693:     int doOutput = 0;
694:     PetscInt bs, b;

696:     if (stateId < 0) {
697:       PetscObjectComposedDataRegister(&stateId);
698:     }
699:     PetscObjectComposedDataGetInt((PetscObject) viewer, stateId, outputState, hasState);
700:     if (!hasState) outputState = 0;
701:     PetscObjectGetName((PetscObject) xin, &name);
702:     VecGetBlockSize(xin, &bs);
703:     if ((bs < 1) || (bs > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %d", bs);
704:     if (format == PETSC_VIEWER_ASCII_VTK) {
705:       if (outputState == 0) {
706:         outputState = 1;
707:         doOutput = 1;
708:       } else if (outputState == 1) doOutput = 0;
709:       else if (outputState == 2) {
710:         outputState = 3;
711:         doOutput = 1;
712:       } else if (outputState == 3) doOutput = 0;
713:       else if (outputState == 4) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");

715:       if (doOutput) {
716:         PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", n/bs);
717:       }
718:     } else {
719:       if (outputState == 0) {
720:         outputState = 2;
721:         doOutput = 1;
722:       } else if (outputState == 1) {
723:         outputState = 4;
724:         doOutput = 1;
725:       } else if (outputState == 2) doOutput = 0;
726:       else if (outputState == 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
727:       else if (outputState == 4) doOutput = 0;

729:       if (doOutput) {
730:         PetscViewerASCIIPrintf(viewer, "CELL_DATA %d\n", n);
731:       }
732:     }
733:     PetscObjectComposedDataSetInt((PetscObject) viewer, stateId, outputState);
734:     if (name) {
735:       if (bs == 3) {
736:         PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
737:       } else {
738:         PetscViewerASCIIPrintf(viewer, "SCALARS %s double %d\n", name, bs);
739:       }
740:     } else {
741:       PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", bs);
742:     }
743:     if (bs != 3) {
744:       PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
745:     }
746:     for (i=0; i<n/bs; i++) {
747:       for (b=0; b<bs; b++) {
748:         if (b > 0) {
749:           PetscViewerASCIIPrintf(viewer," ");
750:         }
751: #if !defined(PETSC_USE_COMPLEX)
752:         PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
753: #endif
754:       }
755:       PetscViewerASCIIPrintf(viewer,"\n");
756:     }
757:   } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS) {
758:     PetscInt bs, b;

760:     VecGetBlockSize(xin, &bs);
761:     if ((bs < 1) || (bs > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "VTK can only handle 3D objects, but vector dimension is %d", bs);
762:     for (i=0; i<n/bs; i++) {
763:       for (b=0; b<bs; b++) {
764:         if (b > 0) {
765:           PetscViewerASCIIPrintf(viewer," ");
766:         }
767: #if !defined(PETSC_USE_COMPLEX)
768:         PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
769: #endif
770:       }
771:       for (b=bs; b<3; b++) {
772:         PetscViewerASCIIPrintf(viewer," 0.0");
773:       }
774:       PetscViewerASCIIPrintf(viewer,"\n");
775:     }
776:   } else if (format == PETSC_VIEWER_ASCII_PCICE) {
777:     PetscInt bs, b;

779:     VecGetBlockSize(xin, &bs);
780:     if ((bs < 1) || (bs > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "PCICE can only handle up to 3D objects, but vector dimension is %d", bs);
781:     PetscViewerASCIIPrintf(viewer,"%D\n", xin->map->N/bs);
782:     for (i=0; i<n/bs; i++) {
783:       PetscViewerASCIIPrintf(viewer,"%7D   ", i+1);
784:       for (b=0; b<bs; b++) {
785:         if (b > 0) {
786:           PetscViewerASCIIPrintf(viewer," ");
787:         }
788: #if !defined(PETSC_USE_COMPLEX)
789:         PetscViewerASCIIPrintf(viewer,"% 12.5E",(double)xv[i*bs+b]);
790: #endif
791:       }
792:       PetscViewerASCIIPrintf(viewer,"\n");
793:     }
794:   } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
795:     /* No info */
796:   } else {
797:     for (i=0; i<n; i++) {
798:       if (format == PETSC_VIEWER_ASCII_INDEX) {
799:         PetscViewerASCIIPrintf(viewer,"%D: ",i);
800:       }
801: #if defined(PETSC_USE_COMPLEX)
802:       if (PetscImaginaryPart(xv[i]) > 0.0) {
803:         PetscViewerASCIIPrintf(viewer,"%g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
804:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
805:         PetscViewerASCIIPrintf(viewer,"%g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
806:       } else {
807:         PetscViewerASCIIPrintf(viewer,"%g\n",(double)PetscRealPart(xv[i]));
808:       }
809: #else
810:       PetscViewerASCIIPrintf(viewer,"%g\n",(double)xv[i]);
811: #endif
812:     }
813:   }
814:   PetscViewerFlush(viewer);
815:   VecRestoreArrayRead(xin,&xv);
816:   return(0);
817: }

819: #include <petscdraw.h>
822: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
823: {
824:   PetscErrorCode    ierr;
825:   PetscInt          i,c,bs = PetscAbs(xin->map->bs),n = xin->map->n/bs;
826:   PetscDraw         win;
827:   PetscReal         *xx;
828:   PetscDrawLG       lg;
829:   const PetscScalar *xv;
830:   PetscReal         *yy;

833:   PetscMalloc1(n,&xx);
834:   PetscMalloc1(n,&yy);
835:   VecGetArrayRead(xin,&xv);
836:   for (c=0; c<bs; c++) {
837:     PetscViewerDrawGetDrawLG(v,c,&lg);
838:     PetscDrawLGGetDraw(lg,&win);
839:     PetscDrawCheckResizedWindow(win);
840:     PetscDrawLGReset(lg);
841:     for (i=0; i<n; i++) {
842:       xx[i] = (PetscReal) i;
843:       yy[i] = PetscRealPart(xv[c + i*bs]);
844:     }
845:     PetscDrawLGAddPoints(lg,n,&xx,&yy);
846:     PetscDrawLGDraw(lg);
847:     PetscDrawSynchronizedFlush(win);
848:   }
849:   VecRestoreArrayRead(xin,&xv);
850:   PetscFree(yy);
851:   PetscFree(xx);
852:   return(0);
853: }

857: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
858: {
859:   PetscErrorCode    ierr;
860:   PetscDraw         draw;
861:   PetscBool         isnull;
862:   PetscViewerFormat format;

865:   PetscViewerDrawGetDraw(v,0,&draw);
866:   PetscDrawIsNull(draw,&isnull); if (isnull) return(0);

868:   PetscViewerGetFormat(v,&format);
869:   /*
870:      Currently it only supports drawing to a line graph */
871:   if (format != PETSC_VIEWER_DRAW_LG) {
872:     PetscViewerPushFormat(v,PETSC_VIEWER_DRAW_LG);
873:   }
874:   VecView_Seq_Draw_LG(xin,v);
875:   if (format != PETSC_VIEWER_DRAW_LG) {
876:     PetscViewerPopFormat(v);
877:   }
878:   return(0);
879: }

883: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
884: {
885:   PetscErrorCode    ierr;
886:   int               fdes;
887:   PetscInt          n = xin->map->n,classid=VEC_FILE_CLASSID;
888:   FILE              *file;
889:   const PetscScalar *xv;
890: #if defined(PETSC_HAVE_MPIIO)
891:   PetscBool         isMPIIO;
892: #endif
893:   PetscBool         skipHeader;
894:   PetscViewerFormat format;

897:   /* Write vector header */
898:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
899:   if (!skipHeader) {
900:     PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
901:     PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
902:   }

904:   /* Write vector contents */
905: #if defined(PETSC_HAVE_MPIIO)
906:   PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);
907:   if (!isMPIIO) {
908: #endif
909:     PetscViewerBinaryGetDescriptor(viewer,&fdes);
910:     VecGetArrayRead(xin,&xv);
911:     PetscBinaryWrite(fdes,(void*)xv,n,PETSC_SCALAR,PETSC_FALSE);
912:     VecRestoreArrayRead(xin,&xv);
913:     PetscViewerGetFormat(viewer,&format);
914:     if (format == PETSC_VIEWER_BINARY_MATLAB) {
915:       MPI_Comm   comm;
916:       FILE       *info;
917:       const char *name;

919:       PetscObjectGetName((PetscObject)xin,&name);
920:       PetscObjectGetComm((PetscObject)viewer,&comm);
921:       PetscViewerBinaryGetInfoPointer(viewer,&info);
922:       PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
923:       PetscFPrintf(comm,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
924:       PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
925:     }
926: #if defined(PETSC_HAVE_MPIIO)
927:   } else {
928:     MPI_Offset   off;
929:     MPI_File     mfdes;
930:     PetscMPIInt  lsize;

932:     PetscMPIIntCast(n,&lsize);
933:     PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
934:     PetscViewerBinaryGetMPIIOOffset(viewer,&off);
935:     MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
936:     VecGetArrayRead(xin,&xv);
937:     MPIU_File_write_all(mfdes,(void*)xv,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
938:     VecRestoreArrayRead(xin,&xv);
939:     PetscViewerBinaryAddMPIIOOffset(viewer,n*sizeof(PetscScalar));
940:   }
941: #endif

943:   PetscViewerBinaryGetInfoPointer(viewer,&file);
944:   if (file) {
945:     if (((PetscObject)xin)->prefix) {
946:       PetscFPrintf(PETSC_COMM_SELF,file,"-%s_vecload_block_size %D\n",((PetscObject)xin)->prefix,PetscAbs(xin->map->bs));
947:     } else {
948:       PetscFPrintf(PETSC_COMM_SELF,file,"-vecload_block_size %D\n",PetscAbs(xin->map->bs));
949:     }
950:   }
951:   return(0);
952: }

954: #if defined(PETSC_HAVE_MATLAB_ENGINE)
955: #include <petscmatlab.h>
956: #include <mat.h>   /* MATLAB include file */
959: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
960: {
961:   PetscErrorCode    ierr;
962:   PetscInt          n;
963:   const PetscScalar *array;

966:   VecGetLocalSize(vec,&n);
967:   PetscObjectName((PetscObject)vec);
968:   VecGetArrayRead(vec,&array);
969:   PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
970:   VecRestoreArrayRead(vec,&array);
971:   return(0);
972: }
973: #endif

977: PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
978: {
980:   PetscBool      isdraw,iascii,issocket,isbinary;
981: #if defined(PETSC_HAVE_MATHEMATICA)
982:   PetscBool      ismathematica;
983: #endif
984: #if defined(PETSC_HAVE_MATLAB_ENGINE)
985:   PetscBool      ismatlab;
986: #endif
987: #if defined(PETSC_HAVE_HDF5)
988:   PetscBool      ishdf5;
989: #endif

992:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
993:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
994:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
995:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
996: #if defined(PETSC_HAVE_MATHEMATICA)
997:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
998: #endif
999: #if defined(PETSC_HAVE_HDF5)
1000:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
1001: #endif
1002: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1003:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
1004: #endif

1006:   if (isdraw) {
1007:     VecView_Seq_Draw(xin,viewer);
1008:   } else if (iascii) {
1009:     VecView_Seq_ASCII(xin,viewer);
1010:   } else if (isbinary) {
1011:     VecView_Seq_Binary(xin,viewer);
1012: #if defined(PETSC_HAVE_MATHEMATICA)
1013:   } else if (ismathematica) {
1014:     PetscViewerMathematicaPutVector(viewer,xin);
1015: #endif
1016: #if defined(PETSC_HAVE_HDF5)
1017:   } else if (ishdf5) {
1018:     VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
1019: #endif
1020: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1021:   } else if (ismatlab) {
1022:     VecView_Seq_Matlab(xin,viewer);
1023: #endif
1024:   }
1025:   return(0);
1026: }

1030: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
1031: {
1032:   const PetscScalar *xx;
1033:   PetscInt          i;
1034:   PetscErrorCode    ierr;

1037:   VecGetArrayRead(xin,&xx);
1038:   for (i=0; i<ni; i++) {
1039:     if (xin->stash.ignorenegidx && ix[i] < 0) continue;
1040: #if defined(PETSC_USE_DEBUG)
1041:     if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
1042:     if (ix[i] >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D to large maximum allowed %D",ix[i],xin->map->n);
1043: #endif
1044:     y[i] = xx[ix[i]];
1045:   }
1046:   VecRestoreArrayRead(xin,&xx);
1047:   return(0);
1048: }

1052: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
1053: {
1054:   PetscScalar    *xx;
1055:   PetscInt       i;

1059:   VecGetArray(xin,&xx);
1060:   if (m == INSERT_VALUES) {
1061:     for (i=0; i<ni; i++) {
1062:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
1063: #if defined(PETSC_USE_DEBUG)
1064:       if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
1065:       if (ix[i] >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",ix[i],xin->map->n);
1066: #endif
1067:       xx[ix[i]] = y[i];
1068:     }
1069:   } else {
1070:     for (i=0; i<ni; i++) {
1071:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
1072: #if defined(PETSC_USE_DEBUG)
1073:       if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
1074:       if (ix[i] >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",ix[i],xin->map->n);
1075: #endif
1076:       xx[ix[i]] += y[i];
1077:     }
1078:   }
1079:   VecRestoreArray(xin,&xx);
1080:   return(0);
1081: }

1085: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
1086: {
1087:   PetscScalar    *xx,*y = (PetscScalar*)yin;
1088:   PetscInt       i,bs,start,j;

1091:   /*
1092:        For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
1093:   */
1095:   VecGetBlockSize(xin,&bs);
1096:   VecGetArray(xin,&xx);
1097:   if (m == INSERT_VALUES) {
1098:     for (i=0; i<ni; i++) {
1099:       start = bs*ix[i];
1100:       if (start < 0) continue;
1101: #if defined(PETSC_USE_DEBUG)
1102:       if (start >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",start,xin->map->n);
1103: #endif
1104:       for (j=0; j<bs; j++) xx[start+j] = y[j];
1105:       y += bs;
1106:     }
1107:   } else {
1108:     for (i=0; i<ni; i++) {
1109:       start = bs*ix[i];
1110:       if (start < 0) continue;
1111: #if defined(PETSC_USE_DEBUG)
1112:       if (start >= xin->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",start,xin->map->n);
1113: #endif
1114:       for (j=0; j<bs; j++) xx[start+j] += y[j];
1115:       y += bs;
1116:     }
1117:   }
1118:   VecRestoreArray(xin,&xx);
1119:   return(0);
1120: }

1124: PetscErrorCode VecDestroy_Seq(Vec v)
1125: {
1126:   Vec_Seq        *vs = (Vec_Seq*)v->data;

1130: #if defined(PETSC_USE_LOG)
1131:   PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
1132: #endif
1133:   PetscFree(vs->array_allocated);
1134:   PetscFree(v->data);
1135:   return(0);
1136: }

1140: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
1141: {
1143:   if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
1144:   return(0);
1145: }

1149: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
1150: {

1154:   VecCreate(PetscObjectComm((PetscObject)win),V);
1155:   PetscObjectSetPrecision((PetscObject)*V,((PetscObject)win)->precision);
1156:   VecSetSizes(*V,win->map->n,win->map->n);
1157:   VecSetType(*V,((PetscObject)win)->type_name);
1158:   PetscLayoutReference(win->map,&(*V)->map);
1159:   PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
1160:   PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);

1162:   (*V)->ops->view          = win->ops->view;
1163:   (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
1164:   return(0);
1165: }

1167: static struct _VecOps DvOps = {VecDuplicate_Seq, /* 1 */
1168:                                VecDuplicateVecs_Default,
1169:                                VecDestroyVecs_Default,
1170:                                VecDot_Seq,
1171:                                VecMDot_Seq,
1172:                                VecNorm_Seq,
1173:                                VecTDot_Seq,
1174:                                VecMTDot_Seq,
1175:                                VecScale_Seq,
1176:                                VecCopy_Seq, /* 10 */
1177:                                VecSet_Seq,
1178:                                VecSwap_Seq,
1179:                                VecAXPY_Seq,
1180:                                VecAXPBY_Seq,
1181:                                VecMAXPY_Seq,
1182:                                VecAYPX_Seq,
1183:                                VecWAXPY_Seq,
1184:                                VecAXPBYPCZ_Seq,
1185:                                VecPointwiseMult_Seq,
1186:                                VecPointwiseDivide_Seq,
1187:                                VecSetValues_Seq, /* 20 */
1188:                                0,0,
1189:                                0,
1190:                                VecGetSize_Seq,
1191:                                VecGetSize_Seq,
1192:                                0,
1193:                                VecMax_Seq,
1194:                                VecMin_Seq,
1195:                                VecSetRandom_Seq,
1196:                                VecSetOption_Seq, /* 30 */
1197:                                VecSetValuesBlocked_Seq,
1198:                                VecDestroy_Seq,
1199:                                VecView_Seq,
1200:                                VecPlaceArray_Seq,
1201:                                VecReplaceArray_Seq,
1202:                                VecDot_Seq,
1203:                                VecTDot_Seq,
1204:                                VecNorm_Seq,
1205:                                VecMDot_Seq,
1206:                                VecMTDot_Seq, /* 40 */
1207:                                VecLoad_Default,
1208:                                VecReciprocal_Default,
1209:                                VecConjugate_Seq,
1210:                                0,
1211:                                0,
1212:                                VecResetArray_Seq,
1213:                                0,
1214:                                VecMaxPointwiseDivide_Seq,
1215:                                VecPointwiseMax_Seq,
1216:                                VecPointwiseMaxAbs_Seq,
1217:                                VecPointwiseMin_Seq,
1218:                                VecGetValues_Seq,
1219:                                0,
1220:                                0,
1221:                                0,
1222:                                0,
1223:                                0,
1224:                                0,
1225:                                VecStrideGather_Default,
1226:                                VecStrideScatter_Default,
1227:                                0,
1228:                                0,
1229:                                0,
1230:                                0,
1231:                                0,
1232:                                VecStrideSubSetGather_Default,
1233:                                VecStrideSubSetScatter_Default
1234: };


1237: /*
1238:       This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
1239: */
1242: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
1243: {
1244:   Vec_Seq        *s;

1248:   PetscNewLog(v,&s);
1249:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));

1251:   v->data            = (void*)s;
1252:   v->petscnative     = PETSC_TRUE;
1253:   s->array           = (PetscScalar*)array;
1254:   s->array_allocated = 0;

1256:   PetscLayoutSetUp(v->map);
1257:   PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
1258: #if defined(PETSC_HAVE_MATLAB_ENGINE)
1259:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
1260:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
1261: #endif
1262: #if defined(PETSC_USE_MIXED_PRECISION)
1263:   ((PetscObject)v)->precision = (PetscPrecision)sizeof(PetscReal);
1264: #endif
1265:   return(0);
1266: }

1270: /*@C
1271:    VecCreateSeqWithArray - Creates a standard,sequential array-style vector,
1272:    where the user provides the array space to store the vector values.

1274:    Collective on MPI_Comm

1276:    Input Parameter:
1277: +  comm - the communicator, should be PETSC_COMM_SELF
1278: .  bs - the block size
1279: .  n - the vector length
1280: -  array - memory where the vector elements are to be stored.

1282:    Output Parameter:
1283: .  V - the vector

1285:    Notes:
1286:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
1287:    same type as an existing vector.

1289:    If the user-provided array is NULL, then VecPlaceArray() can be used
1290:    at a later stage to SET the array for storing the vector values.

1292:    PETSc does NOT free the array when the vector is destroyed via VecDestroy().
1293:    The user should not free the array until the vector is destroyed.

1295:    Level: intermediate

1297:    Concepts: vectors^creating with array

1299: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
1300:           VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
1301: @*/
1302: PetscErrorCode  VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
1303: {
1305:   PetscMPIInt    size;

1308:   VecCreate(comm,V);
1309:   VecSetSizes(*V,n,n);
1310:   VecSetBlockSize(*V,bs);
1311:   MPI_Comm_size(comm,&size);
1312:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
1313:   VecCreate_Seq_Private(*V,array);
1314:   return(0);
1315: }