Actual source code: bvec2.c

petsc-3.7.3 2016-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>

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

 16: PetscErrorCode VecPointwiseMax_Seq(Vec win,Vec xin,Vec yin)
 17: {
 19:   PetscInt       n = win->map->n,i;
 20:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

 23:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
 24:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
 25:   VecGetArray(win,&ww);

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

 29:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
 30:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
 31:   VecRestoreArray(win,&ww);
 32:   PetscLogFlops(n);
 33:   return(0);
 34: }

 38: PetscErrorCode VecPointwiseMin_Seq(Vec win,Vec xin,Vec yin)
 39: {
 41:   PetscInt       n = win->map->n,i;
 42:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

 45:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
 46:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
 47:   VecGetArray(win,&ww);

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

 51:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
 52:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
 53:   VecRestoreArray(win,&ww);
 54:   PetscLogFlops(n);
 55:   return(0);
 56: }

 60: PetscErrorCode VecPointwiseMaxAbs_Seq(Vec win,Vec xin,Vec yin)
 61: {
 63:   PetscInt       n = win->map->n,i;
 64:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

 67:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
 68:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
 69:   VecGetArray(win,&ww);

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

 73:   PetscLogFlops(n);
 74:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
 75:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
 76:   VecRestoreArray(win,&ww);
 77:   return(0);
 78: }

 80: #include <../src/vec/vec/impls/seq/ftn-kernels/fxtimesy.h>

 84: PetscErrorCode VecPointwiseMult_Seq(Vec win,Vec xin,Vec yin)
 85: {
 87:   PetscInt       n = win->map->n,i;
 88:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

 91:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
 92:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
 93:   VecGetArray(win,&ww);
 94:   if (ww == xx) {
 95:     for (i=0; i<n; i++) ww[i] *= yy[i];
 96:   } else if (ww == yy) {
 97:     for (i=0; i<n; i++) ww[i] *= xx[i];
 98:   } else {
 99: #if defined(PETSC_USE_FORTRAN_KERNEL_XTIMESY)
100:     fortranxtimesy_(xx,yy,ww,&n);
101: #else
102:     for (i=0; i<n; i++) ww[i] = xx[i] * yy[i];
103: #endif
104:   }
105:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
106:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
107:   VecRestoreArray(win,&ww);
108:   PetscLogFlops(n);
109:   return(0);
110: }

114: PetscErrorCode VecPointwiseDivide_Seq(Vec win,Vec xin,Vec yin)
115: {
117:   PetscInt       n = win->map->n,i;
118:   PetscScalar    *ww,*xx,*yy; /* cannot make xx or yy const since might be ww */

121:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
122:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
123:   VecGetArray(win,&ww);

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

127:   PetscLogFlops(n);
128:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
129:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
130:   VecRestoreArray(win,&ww);
131:   return(0);
132: }

136: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
137: {
139:   PetscInt       n = xin->map->n,i;
140:   PetscScalar    *xx;

143:   VecGetArray(xin,&xx);
144:   for (i=0; i<n; i++) {PetscRandomGetValue(r,&xx[i]);}
145:   VecRestoreArray(xin,&xx);
146:   return(0);
147: }

151: PetscErrorCode VecGetSize_Seq(Vec vin,PetscInt *size)
152: {
154:   *size = vin->map->n;
155:   return(0);
156: }



162: PetscErrorCode VecConjugate_Seq(Vec xin)
163: {
164:   PetscScalar    *x;
165:   PetscInt       n = xin->map->n;

169:   VecGetArray(xin,&x);
170:   while (n-->0) {
171:     *x = PetscConj(*x);
172:     x++;
173:   }
174:   VecRestoreArray(xin,&x);
175:   return(0);
176: }

180: PetscErrorCode VecResetArray_Seq(Vec vin)
181: {
182:   Vec_Seq *v = (Vec_Seq*)vin->data;

185:   v->array         = v->unplacedarray;
186:   v->unplacedarray = 0;
187:   return(0);
188: }

192: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
193: {
194:   PetscScalar       *ya;
195:   const PetscScalar *xa;
196:   PetscErrorCode    ierr;

199:   if (xin != yin) {
200:     VecGetArrayRead(xin,&xa);
201:     VecGetArray(yin,&ya);
202:     PetscMemcpy(ya,xa,xin->map->n*sizeof(PetscScalar));
203:     VecRestoreArrayRead(xin,&xa);
204:     VecRestoreArray(yin,&ya);
205:   }
206:   return(0);
207: }

211: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
212: {
213:   PetscScalar    *ya, *xa;
215:   PetscBLASInt   one = 1,bn;

218:   if (xin != yin) {
219:     PetscBLASIntCast(xin->map->n,&bn);
220:     VecGetArray(xin,&xa);
221:     VecGetArray(yin,&ya);
222:     PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa,&one,ya,&one));
223:     VecRestoreArray(xin,&xa);
224:     VecRestoreArray(yin,&ya);
225:   }
226:   return(0);
227: }

229: #include <../src/vec/vec/impls/seq/ftn-kernels/fnorm.h>

233: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
234: {
235:   const PetscScalar *xx;
236:   PetscErrorCode    ierr;
237:   PetscInt          n = xin->map->n;
238:   PetscBLASInt      one = 1, bn;

241:   PetscBLASIntCast(n,&bn);
242:   if (type == NORM_2 || type == NORM_FROBENIUS) {
243:     VecGetArrayRead(xin,&xx);
244:     *z   = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
245:     *z   = PetscSqrtReal(*z);
246:     VecRestoreArrayRead(xin,&xx);
247:     PetscLogFlops(PetscMax(2.0*n-1,0.0));
248:   } else if (type == NORM_INFINITY) {
249:     PetscInt  i;
250:     PetscReal max = 0.0,tmp;

252:     VecGetArrayRead(xin,&xx);
253:     for (i=0; i<n; i++) {
254:       if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
255:       /* check special case of tmp == NaN */
256:       if (tmp != tmp) {max = tmp; break;}
257:       xx++;
258:     }
259:     VecRestoreArrayRead(xin,&xx);
260:     *z   = max;
261:   } else if (type == NORM_1) {
262: #if defined(PETSC_USE_COMPLEX)
263:     PetscReal tmp = 0.0;
264:     PetscInt    i;
265: #endif
266:     VecGetArrayRead(xin,&xx);
267: #if defined(PETSC_USE_COMPLEX)
268:     /* BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we provide a custom loop instead */
269:     for (i=0; i<n; i++) {
270:       tmp += PetscAbsScalar(xx[i]);
271:     }
272:     *z = tmp;
273: #else
274:     PetscStackCallBLAS("BLASasum",*z   = BLASasum_(&bn,xx,&one));
275: #endif
276:     VecRestoreArrayRead(xin,&xx);
277:     PetscLogFlops(PetscMax(n-1.0,0.0));
278:   } else if (type == NORM_1_AND_2) {
279:     VecNorm_Seq(xin,NORM_1,z);
280:     VecNorm_Seq(xin,NORM_2,z+1);
281:   }
282:   return(0);
283: }

287: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
288: {
289:   PetscErrorCode    ierr;
290:   PetscInt          i,n = xin->map->n;
291:   const char        *name;
292:   PetscViewerFormat format;
293:   const PetscScalar *xv;

296:   VecGetArrayRead(xin,&xv);
297:   PetscViewerGetFormat(viewer,&format);
298:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
299:     PetscObjectGetName((PetscObject)xin,&name);
300:     PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
301:     for (i=0; i<n; i++) {
302: #if defined(PETSC_USE_COMPLEX)
303:       if (PetscImaginaryPart(xv[i]) > 0.0) {
304:         PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
305:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
306:         PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
307:       } else {
308:         PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)PetscRealPart(xv[i]));
309:       }
310: #else
311:       PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double) xv[i]);
312: #endif
313:     }
314:     PetscViewerASCIIPrintf(viewer,"];\n");
315:   } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
316:     for (i=0; i<n; i++) {
317: #if defined(PETSC_USE_COMPLEX)
318:       PetscViewerASCIIPrintf(viewer,"%18.16e %18.16e\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
319: #else
320:       PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)xv[i]);
321: #endif
322:     }
323:   } else if (format == PETSC_VIEWER_ASCII_VTK || format == PETSC_VIEWER_ASCII_VTK_CELL) {
324:     /*
325:        state 0: No header has been output
326:        state 1: Only POINT_DATA has been output
327:        state 2: Only CELL_DATA has been output
328:        state 3: Output both, POINT_DATA last
329:        state 4: Output both, CELL_DATA last
330:     */
331:     static PetscInt stateId = -1;
332:     int outputState = 0;
333:     PetscBool  hasState;
334:     int doOutput = 0;
335:     PetscInt bs, b;

337:     if (stateId < 0) {
338:       PetscObjectComposedDataRegister(&stateId);
339:     }
340:     PetscObjectComposedDataGetInt((PetscObject) viewer, stateId, outputState, hasState);
341:     if (!hasState) outputState = 0;
342:     PetscObjectGetName((PetscObject) xin, &name);
343:     VecGetBlockSize(xin, &bs);
344:     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);
345:     if (format == PETSC_VIEWER_ASCII_VTK) {
346:       if (outputState == 0) {
347:         outputState = 1;
348:         doOutput = 1;
349:       } else if (outputState == 1) doOutput = 0;
350:       else if (outputState == 2) {
351:         outputState = 3;
352:         doOutput = 1;
353:       } else if (outputState == 3) doOutput = 0;
354:       else if (outputState == 4) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");

356:       if (doOutput) {
357:         PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", n/bs);
358:       }
359:     } else {
360:       if (outputState == 0) {
361:         outputState = 2;
362:         doOutput = 1;
363:       } else if (outputState == 1) {
364:         outputState = 4;
365:         doOutput = 1;
366:       } else if (outputState == 2) doOutput = 0;
367:       else if (outputState == 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
368:       else if (outputState == 4) doOutput = 0;

370:       if (doOutput) {
371:         PetscViewerASCIIPrintf(viewer, "CELL_DATA %d\n", n);
372:       }
373:     }
374:     PetscObjectComposedDataSetInt((PetscObject) viewer, stateId, outputState);
375:     if (name) {
376:       if (bs == 3) {
377:         PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
378:       } else {
379:         PetscViewerASCIIPrintf(viewer, "SCALARS %s double %d\n", name, bs);
380:       }
381:     } else {
382:       PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", bs);
383:     }
384:     if (bs != 3) {
385:       PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
386:     }
387:     for (i=0; i<n/bs; i++) {
388:       for (b=0; b<bs; b++) {
389:         if (b > 0) {
390:           PetscViewerASCIIPrintf(viewer," ");
391:         }
392: #if !defined(PETSC_USE_COMPLEX)
393:         PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
394: #endif
395:       }
396:       PetscViewerASCIIPrintf(viewer,"\n");
397:     }
398:   } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS) {
399:     PetscInt bs, b;

401:     VecGetBlockSize(xin, &bs);
402:     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);
403:     for (i=0; i<n/bs; i++) {
404:       for (b=0; b<bs; b++) {
405:         if (b > 0) {
406:           PetscViewerASCIIPrintf(viewer," ");
407:         }
408: #if !defined(PETSC_USE_COMPLEX)
409:         PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
410: #endif
411:       }
412:       for (b=bs; b<3; b++) {
413:         PetscViewerASCIIPrintf(viewer," 0.0");
414:       }
415:       PetscViewerASCIIPrintf(viewer,"\n");
416:     }
417:   } else if (format == PETSC_VIEWER_ASCII_PCICE) {
418:     PetscInt bs, b;

420:     VecGetBlockSize(xin, &bs);
421:     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);
422:     PetscViewerASCIIPrintf(viewer,"%D\n", xin->map->N/bs);
423:     for (i=0; i<n/bs; i++) {
424:       PetscViewerASCIIPrintf(viewer,"%7D   ", i+1);
425:       for (b=0; b<bs; b++) {
426:         if (b > 0) {
427:           PetscViewerASCIIPrintf(viewer," ");
428:         }
429: #if !defined(PETSC_USE_COMPLEX)
430:         PetscViewerASCIIPrintf(viewer,"% 12.5E",(double)xv[i*bs+b]);
431: #endif
432:       }
433:       PetscViewerASCIIPrintf(viewer,"\n");
434:     }
435:   } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
436:     /* No info */
437:   } else {
438:     for (i=0; i<n; i++) {
439:       if (format == PETSC_VIEWER_ASCII_INDEX) {
440:         PetscViewerASCIIPrintf(viewer,"%D: ",i);
441:       }
442: #if defined(PETSC_USE_COMPLEX)
443:       if (PetscImaginaryPart(xv[i]) > 0.0) {
444:         PetscViewerASCIIPrintf(viewer,"%g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
445:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
446:         PetscViewerASCIIPrintf(viewer,"%g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
447:       } else {
448:         PetscViewerASCIIPrintf(viewer,"%g\n",(double)PetscRealPart(xv[i]));
449:       }
450: #else
451:       PetscViewerASCIIPrintf(viewer,"%g\n",(double)xv[i]);
452: #endif
453:     }
454:   }
455:   PetscViewerFlush(viewer);
456:   VecRestoreArrayRead(xin,&xv);
457:   return(0);
458: }

460: #include <petscdraw.h>
463: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
464: {
465:   PetscDraw         draw;
466:   PetscBool         isnull;
467:   PetscDrawLG       lg;
468:   PetscErrorCode    ierr;
469:   PetscInt          i,c,bs = PetscAbs(xin->map->bs),n = xin->map->n/bs;
470:   const PetscScalar *xv;
471:   PetscReal         *xx,*yy;
472:   int               colors[] = {PETSC_DRAW_RED};

475:   PetscViewerDrawGetDraw(v,0,&draw);
476:   PetscDrawIsNull(draw,&isnull);
477:   if (isnull) return(0);

479:   PetscMalloc2(n,&xx,n,&yy);
480:   VecGetArrayRead(xin,&xv);
481:   for (c=0; c<bs; c++) {
482:     PetscViewerDrawGetDrawLG(v,c,&lg);
483:     PetscDrawLGReset(lg);
484:     PetscDrawLGSetDimension(lg,1);
485:     PetscDrawLGSetColors(lg,colors);
486:     for (i=0; i<n; i++) {
487:       xx[i] = (PetscReal)i;
488:       yy[i] = PetscRealPart(xv[c + i*bs]);
489:     }
490:     PetscDrawLGAddPoints(lg,n,&xx,&yy);
491:     PetscDrawLGDraw(lg);
492:     PetscDrawLGSave(lg);
493:   }
494:   VecRestoreArrayRead(xin,&xv);
495:   PetscFree2(xx,yy);
496:   return(0);
497: }

501: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
502: {
503:   PetscErrorCode    ierr;
504:   PetscDraw         draw;
505:   PetscBool         isnull;

508:   PetscViewerDrawGetDraw(v,0,&draw);
509:   PetscDrawIsNull(draw,&isnull);
510:   if (isnull) return(0);
511:   PetscViewerPushFormat(v,PETSC_VIEWER_DRAW_LG);
512:   VecView_Seq_Draw_LG(xin,v);
513:   PetscViewerPopFormat(v);
514:   return(0);
515: }

519: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
520: {
521:   PetscErrorCode    ierr;
522:   int               fdes;
523:   PetscInt          n = xin->map->n,classid=VEC_FILE_CLASSID;
524:   FILE              *file;
525:   const PetscScalar *xv;
526: #if defined(PETSC_HAVE_MPIIO)
527:   PetscBool         isMPIIO;
528: #endif
529:   PetscBool         skipHeader;
530:   PetscViewerFormat format;

533:   /* Write vector header */
534:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
535:   if (!skipHeader) {
536:     PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
537:     PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
538:   }

540:   /* Write vector contents */
541: #if defined(PETSC_HAVE_MPIIO)
542:   PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);
543:   if (!isMPIIO) {
544: #endif
545:     PetscViewerBinaryGetDescriptor(viewer,&fdes);
546:     VecGetArrayRead(xin,&xv);
547:     PetscBinaryWrite(fdes,(void*)xv,n,PETSC_SCALAR,PETSC_FALSE);
548:     VecRestoreArrayRead(xin,&xv);
549:     PetscViewerGetFormat(viewer,&format);
550:     if (format == PETSC_VIEWER_BINARY_MATLAB) {
551:       MPI_Comm   comm;
552:       FILE       *info;
553:       const char *name;

555:       PetscObjectGetName((PetscObject)xin,&name);
556:       PetscObjectGetComm((PetscObject)viewer,&comm);
557:       PetscViewerBinaryGetInfoPointer(viewer,&info);
558:       PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
559:       PetscFPrintf(comm,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
560:       PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
561:     }
562: #if defined(PETSC_HAVE_MPIIO)
563:   } else {
564:     MPI_Offset   off;
565:     MPI_File     mfdes;
566:     PetscMPIInt  lsize;

568:     PetscMPIIntCast(n,&lsize);
569:     PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
570:     PetscViewerBinaryGetMPIIOOffset(viewer,&off);
571:     MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
572:     VecGetArrayRead(xin,&xv);
573:     MPIU_File_write_all(mfdes,(void*)xv,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
574:     VecRestoreArrayRead(xin,&xv);
575:     PetscViewerBinaryAddMPIIOOffset(viewer,n*sizeof(PetscScalar));
576:   }
577: #endif

579:   PetscViewerBinaryGetInfoPointer(viewer,&file);
580:   if (file) {
581:     if (((PetscObject)xin)->prefix) {
582:       PetscFPrintf(PETSC_COMM_SELF,file,"-%svecload_block_size %D\n",((PetscObject)xin)->prefix,PetscAbs(xin->map->bs));
583:     } else {
584:       PetscFPrintf(PETSC_COMM_SELF,file,"-vecload_block_size %D\n",PetscAbs(xin->map->bs));
585:     }
586:   }
587:   return(0);
588: }

590: #if defined(PETSC_HAVE_MATLAB_ENGINE)
591: #include <petscmatlab.h>
592: #include <mat.h>   /* MATLAB include file */
595: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
596: {
597:   PetscErrorCode    ierr;
598:   PetscInt          n;
599:   const PetscScalar *array;

602:   VecGetLocalSize(vec,&n);
603:   PetscObjectName((PetscObject)vec);
604:   VecGetArrayRead(vec,&array);
605:   PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
606:   VecRestoreArrayRead(vec,&array);
607:   return(0);
608: }
609: #endif

613: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
614: {
616:   PetscBool      isdraw,iascii,issocket,isbinary;
617: #if defined(PETSC_HAVE_MATHEMATICA)
618:   PetscBool      ismathematica;
619: #endif
620: #if defined(PETSC_HAVE_MATLAB_ENGINE)
621:   PetscBool      ismatlab;
622: #endif
623: #if defined(PETSC_HAVE_HDF5)
624:   PetscBool      ishdf5;
625: #endif

628:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
629:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
630:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
631:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
632: #if defined(PETSC_HAVE_MATHEMATICA)
633:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
634: #endif
635: #if defined(PETSC_HAVE_HDF5)
636:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
637: #endif
638: #if defined(PETSC_HAVE_MATLAB_ENGINE)
639:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
640: #endif

642:   if (isdraw) {
643:     VecView_Seq_Draw(xin,viewer);
644:   } else if (iascii) {
645:     VecView_Seq_ASCII(xin,viewer);
646:   } else if (isbinary) {
647:     VecView_Seq_Binary(xin,viewer);
648: #if defined(PETSC_HAVE_MATHEMATICA)
649:   } else if (ismathematica) {
650:     PetscViewerMathematicaPutVector(viewer,xin);
651: #endif
652: #if defined(PETSC_HAVE_HDF5)
653:   } else if (ishdf5) {
654:     VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
655: #endif
656: #if defined(PETSC_HAVE_MATLAB_ENGINE)
657:   } else if (ismatlab) {
658:     VecView_Seq_Matlab(xin,viewer);
659: #endif
660:   }
661:   return(0);
662: }

666: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
667: {
668:   const PetscScalar *xx;
669:   PetscInt          i;
670:   PetscErrorCode    ierr;

673:   VecGetArrayRead(xin,&xx);
674:   for (i=0; i<ni; i++) {
675:     if (xin->stash.ignorenegidx && ix[i] < 0) continue;
676: #if defined(PETSC_USE_DEBUG)
677:     if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
678:     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);
679: #endif
680:     y[i] = xx[ix[i]];
681:   }
682:   VecRestoreArrayRead(xin,&xx);
683:   return(0);
684: }

688: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
689: {
690:   PetscScalar    *xx;
691:   PetscInt       i;

695:   VecGetArray(xin,&xx);
696:   if (m == INSERT_VALUES) {
697:     for (i=0; i<ni; i++) {
698:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
699: #if defined(PETSC_USE_DEBUG)
700:       if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
701:       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);
702: #endif
703:       xx[ix[i]] = y[i];
704:     }
705:   } else {
706:     for (i=0; i<ni; i++) {
707:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
708: #if defined(PETSC_USE_DEBUG)
709:       if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
710:       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);
711: #endif
712:       xx[ix[i]] += y[i];
713:     }
714:   }
715:   VecRestoreArray(xin,&xx);
716:   return(0);
717: }

721: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
722: {
723:   PetscScalar    *xx,*y = (PetscScalar*)yin;
724:   PetscInt       i,bs,start,j;

727:   /*
728:        For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
729:   */
731:   VecGetBlockSize(xin,&bs);
732:   VecGetArray(xin,&xx);
733:   if (m == INSERT_VALUES) {
734:     for (i=0; i<ni; i++) {
735:       start = bs*ix[i];
736:       if (start < 0) continue;
737: #if defined(PETSC_USE_DEBUG)
738:       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);
739: #endif
740:       for (j=0; j<bs; j++) xx[start+j] = y[j];
741:       y += bs;
742:     }
743:   } else {
744:     for (i=0; i<ni; i++) {
745:       start = bs*ix[i];
746:       if (start < 0) continue;
747: #if defined(PETSC_USE_DEBUG)
748:       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);
749: #endif
750:       for (j=0; j<bs; j++) xx[start+j] += y[j];
751:       y += bs;
752:     }
753:   }
754:   VecRestoreArray(xin,&xx);
755:   return(0);
756: }

760: PetscErrorCode VecDestroy_Seq(Vec v)
761: {
762:   Vec_Seq        *vs = (Vec_Seq*)v->data;

766: #if defined(PETSC_USE_LOG)
767:   PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
768: #endif
769:   PetscFree(vs->array_allocated);
770:   PetscFree(v->data);
771:   return(0);
772: }

776: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
777: {
779:   if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
780:   return(0);
781: }

785: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
786: {

790:   VecCreate(PetscObjectComm((PetscObject)win),V);
791:   PetscObjectSetPrecision((PetscObject)*V,((PetscObject)win)->precision);
792:   VecSetSizes(*V,win->map->n,win->map->n);
793:   VecSetType(*V,((PetscObject)win)->type_name);
794:   PetscLayoutReference(win->map,&(*V)->map);
795:   PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
796:   PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);

798:   (*V)->ops->view          = win->ops->view;
799:   (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
800:   return(0);
801: }

803: static struct _VecOps DvOps = {VecDuplicate_Seq, /* 1 */
804:                                VecDuplicateVecs_Default,
805:                                VecDestroyVecs_Default,
806:                                VecDot_Seq,
807:                                VecMDot_Seq,
808:                                VecNorm_Seq,
809:                                VecTDot_Seq,
810:                                VecMTDot_Seq,
811:                                VecScale_Seq,
812:                                VecCopy_Seq, /* 10 */
813:                                VecSet_Seq,
814:                                VecSwap_Seq,
815:                                VecAXPY_Seq,
816:                                VecAXPBY_Seq,
817:                                VecMAXPY_Seq,
818:                                VecAYPX_Seq,
819:                                VecWAXPY_Seq,
820:                                VecAXPBYPCZ_Seq,
821:                                VecPointwiseMult_Seq,
822:                                VecPointwiseDivide_Seq,
823:                                VecSetValues_Seq, /* 20 */
824:                                0,0,
825:                                0,
826:                                VecGetSize_Seq,
827:                                VecGetSize_Seq,
828:                                0,
829:                                VecMax_Seq,
830:                                VecMin_Seq,
831:                                VecSetRandom_Seq,
832:                                VecSetOption_Seq, /* 30 */
833:                                VecSetValuesBlocked_Seq,
834:                                VecDestroy_Seq,
835:                                VecView_Seq,
836:                                VecPlaceArray_Seq,
837:                                VecReplaceArray_Seq,
838:                                VecDot_Seq,
839:                                VecTDot_Seq,
840:                                VecNorm_Seq,
841:                                VecMDot_Seq,
842:                                VecMTDot_Seq, /* 40 */
843:                                VecLoad_Default,
844:                                VecReciprocal_Default,
845:                                VecConjugate_Seq,
846:                                0,
847:                                0,
848:                                VecResetArray_Seq,
849:                                0,
850:                                VecMaxPointwiseDivide_Seq,
851:                                VecPointwiseMax_Seq,
852:                                VecPointwiseMaxAbs_Seq,
853:                                VecPointwiseMin_Seq,
854:                                VecGetValues_Seq,
855:                                0,
856:                                0,
857:                                0,
858:                                0,
859:                                0,
860:                                0,
861:                                VecStrideGather_Default,
862:                                VecStrideScatter_Default,
863:                                0,
864:                                0,
865:                                0,
866:                                0,
867:                                0,
868:                                VecStrideSubSetGather_Default,
869:                                VecStrideSubSetScatter_Default,
870:                                0,
871:                                0
872: };


875: /*
876:       This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
877: */
880: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
881: {
882:   Vec_Seq        *s;

886:   PetscNewLog(v,&s);
887:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));

889:   v->data            = (void*)s;
890:   v->petscnative     = PETSC_TRUE;
891:   s->array           = (PetscScalar*)array;
892:   s->array_allocated = 0;

894:   PetscLayoutSetUp(v->map);
895:   PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
896: #if defined(PETSC_HAVE_MATLAB_ENGINE)
897:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
898:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
899: #endif
900: #if defined(PETSC_USE_MIXED_PRECISION)
901:   ((PetscObject)v)->precision = (PetscPrecision)sizeof(PetscReal);
902: #endif
903:   return(0);
904: }

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

912:    Collective on MPI_Comm

914:    Input Parameter:
915: +  comm - the communicator, should be PETSC_COMM_SELF
916: .  bs - the block size
917: .  n - the vector length
918: -  array - memory where the vector elements are to be stored.

920:    Output Parameter:
921: .  V - the vector

923:    Notes:
924:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
925:    same type as an existing vector.

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

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

933:    Level: intermediate

935:    Concepts: vectors^creating with array

937: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
938:           VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
939: @*/
940: PetscErrorCode  VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
941: {
943:   PetscMPIInt    size;

946:   VecCreate(comm,V);
947:   VecSetSizes(*V,n,n);
948:   VecSetBlockSize(*V,bs);
949:   MPI_Comm_size(comm,&size);
950:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
951:   VecCreate_Seq_Private(*V,array);
952:   return(0);
953: }