Actual source code: bvec2.c

petsc-master 2019-09-17
Report Typos and Errors

  2: /*
  3:    Implements the sequential vectors.
  4: */

  6:  #include <../src/vec/vec/impls/dvecimpl.h>
  7: #include <../src/vec/vec/impls/mpi/pvecimpl.h>      /* For VecView_MPI_HDF5 */
  8:  #include <petsc/private/glvisviewerimpl.h>
  9:  #include <petsc/private/glvisvecimpl.h>
 10:  #include <petscblaslapack.h>

 12: #if defined(PETSC_HAVE_HDF5)
 13: extern PetscErrorCode VecView_MPI_HDF5(Vec,PetscViewer);
 14: #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: }

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

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

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

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

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

 63:   VecGetArrayRead(xin,(const PetscScalar**)&xx);
 64:   VecGetArrayRead(yin,(const PetscScalar**)&yy);
 65:   VecGetArray(win,&ww);

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

 69:   PetscLogFlops(n);
 70:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
 71:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
 72:   VecRestoreArray(win,&ww);
 73:   return(0);
 74: }

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

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

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

106: PetscErrorCode VecPointwiseDivide_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++) {
118:     if (yy[i] != 0.0) ww[i] = xx[i] / yy[i];
119:     else ww[i] = 0.0;
120:   }

122:   PetscLogFlops(n);
123:   VecRestoreArrayRead(xin,(const PetscScalar**)&xx);
124:   VecRestoreArrayRead(yin,(const PetscScalar**)&yy);
125:   VecRestoreArray(win,&ww);
126:   return(0);
127: }

129: PetscErrorCode VecSetRandom_Seq(Vec xin,PetscRandom r)
130: {
132:   PetscInt       n = xin->map->n,i;
133:   PetscScalar    *xx;

136:   VecGetArray(xin,&xx);
137:   for (i=0; i<n; i++) {PetscRandomGetValue(r,&xx[i]);}
138:   VecRestoreArray(xin,&xx);
139:   return(0);
140: }

142: PetscErrorCode VecGetSize_Seq(Vec vin,PetscInt *size)
143: {
145:   *size = vin->map->n;
146:   return(0);
147: }



151: PetscErrorCode VecConjugate_Seq(Vec xin)
152: {
153:   PetscScalar    *x;
154:   PetscInt       n = xin->map->n;

158:   VecGetArray(xin,&x);
159:   while (n-->0) {
160:     *x = PetscConj(*x);
161:     x++;
162:   }
163:   VecRestoreArray(xin,&x);
164:   return(0);
165: }

167: PetscErrorCode VecResetArray_Seq(Vec vin)
168: {
169:   Vec_Seq *v = (Vec_Seq*)vin->data;

172:   v->array         = v->unplacedarray;
173:   v->unplacedarray = 0;
174:   return(0);
175: }

177: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
178: {
179:   PetscScalar       *ya;
180:   const PetscScalar *xa;
181:   PetscErrorCode    ierr;

184:   if (xin != yin) {
185:     VecGetArrayRead(xin,&xa);
186:     VecGetArray(yin,&ya);
187:     PetscArraycpy(ya,xa,xin->map->n);
188:     VecRestoreArrayRead(xin,&xa);
189:     VecRestoreArray(yin,&ya);
190:   }
191:   return(0);
192: }

194: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
195: {
196:   PetscScalar    *ya, *xa;
198:   PetscBLASInt   one = 1,bn;

201:   if (xin != yin) {
202:     PetscBLASIntCast(xin->map->n,&bn);
203:     VecGetArray(xin,&xa);
204:     VecGetArray(yin,&ya);
205:     PetscStackCallBLAS("BLASswap",BLASswap_(&bn,xa,&one,ya,&one));
206:     VecRestoreArray(xin,&xa);
207:     VecRestoreArray(yin,&ya);
208:   }
209:   return(0);
210: }

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

214: PetscErrorCode VecNorm_Seq(Vec xin,NormType type,PetscReal *z)
215: {
216:   const PetscScalar *xx;
217:   PetscErrorCode    ierr;
218:   PetscInt          n = xin->map->n;
219:   PetscBLASInt      one = 1, bn;

222:   PetscBLASIntCast(n,&bn);
223:   if (type == NORM_2 || type == NORM_FROBENIUS) {
224:     VecGetArrayRead(xin,&xx);
225: #if defined(PETSC_USE_REAL___FP16)
226:     *z   = BLASnrm2_(&bn,xx,&one);
227: #else
228:     *z   = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
229:     *z   = PetscSqrtReal(*z);
230: #endif
231:     VecRestoreArrayRead(xin,&xx);
232:     PetscLogFlops(PetscMax(2.0*n-1,0.0));
233:   } else if (type == NORM_INFINITY) {
234:     PetscInt  i;
235:     PetscReal max = 0.0,tmp;

237:     VecGetArrayRead(xin,&xx);
238:     for (i=0; i<n; i++) {
239:       if ((tmp = PetscAbsScalar(*xx)) > max) max = tmp;
240:       /* check special case of tmp == NaN */
241:       if (tmp != tmp) {max = tmp; break;}
242:       xx++;
243:     }
244:     VecRestoreArrayRead(xin,&xx);
245:     *z   = max;
246:   } else if (type == NORM_1) {
247: #if defined(PETSC_USE_COMPLEX)
248:     PetscReal tmp = 0.0;
249:     PetscInt    i;
250: #endif
251:     VecGetArrayRead(xin,&xx);
252: #if defined(PETSC_USE_COMPLEX)
253:     /* BLASasum() returns the nonstandard 1 norm of the 1 norm of the complex entries so we provide a custom loop instead */
254:     for (i=0; i<n; i++) {
255:       tmp += PetscAbsScalar(xx[i]);
256:     }
257:     *z = tmp;
258: #else
259:     PetscStackCallBLAS("BLASasum",*z   = BLASasum_(&bn,xx,&one));
260: #endif
261:     VecRestoreArrayRead(xin,&xx);
262:     PetscLogFlops(PetscMax(n-1.0,0.0));
263:   } else if (type == NORM_1_AND_2) {
264:     VecNorm_Seq(xin,NORM_1,z);
265:     VecNorm_Seq(xin,NORM_2,z+1);
266:   }
267:   return(0);
268: }

270: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
271: {
272:   PetscErrorCode    ierr;
273:   PetscInt          i,n = xin->map->n;
274:   const char        *name;
275:   PetscViewerFormat format;
276:   const PetscScalar *xv;

279:   VecGetArrayRead(xin,&xv);
280:   PetscViewerGetFormat(viewer,&format);
281:   if (format == PETSC_VIEWER_ASCII_MATLAB) {
282:     PetscObjectGetName((PetscObject)xin,&name);
283:     PetscViewerASCIIPrintf(viewer,"%s = [\n",name);
284:     for (i=0; i<n; i++) {
285: #if defined(PETSC_USE_COMPLEX)
286:       if (PetscImaginaryPart(xv[i]) > 0.0) {
287:         PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16ei\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
288:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
289:         PetscViewerASCIIPrintf(viewer,"%18.16e - %18.16ei\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
290:       } else {
291:         PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)PetscRealPart(xv[i]));
292:       }
293: #else
294:       PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double) xv[i]);
295: #endif
296:     }
297:     PetscViewerASCIIPrintf(viewer,"];\n");
298:   } else if (format == PETSC_VIEWER_ASCII_SYMMODU) {
299:     for (i=0; i<n; i++) {
300: #if defined(PETSC_USE_COMPLEX)
301:       PetscViewerASCIIPrintf(viewer,"%18.16e %18.16e\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
302: #else
303:       PetscViewerASCIIPrintf(viewer,"%18.16e\n",(double)xv[i]);
304: #endif
305:     }
306:   } else if (format == PETSC_VIEWER_ASCII_VTK || format == PETSC_VIEWER_ASCII_VTK_CELL) {
307:     /*
308:        state 0: No header has been output
309:        state 1: Only POINT_DATA has been output
310:        state 2: Only CELL_DATA has been output
311:        state 3: Output both, POINT_DATA last
312:        state 4: Output both, CELL_DATA last
313:     */
314:     static PetscInt stateId = -1;
315:     int outputState = 0;
316:     PetscBool  hasState;
317:     int doOutput = 0;
318:     PetscInt bs, b;

320:     if (stateId < 0) {
321:       PetscObjectComposedDataRegister(&stateId);
322:     }
323:     PetscObjectComposedDataGetInt((PetscObject) viewer, stateId, outputState, hasState);
324:     if (!hasState) outputState = 0;
325:     PetscObjectGetName((PetscObject) xin, &name);
326:     VecGetBlockSize(xin, &bs);
327:     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);
328:     if (format == PETSC_VIEWER_ASCII_VTK) {
329:       if (outputState == 0) {
330:         outputState = 1;
331:         doOutput = 1;
332:       } else if (outputState == 1) doOutput = 0;
333:       else if (outputState == 2) {
334:         outputState = 3;
335:         doOutput = 1;
336:       } else if (outputState == 3) doOutput = 0;
337:       else if (outputState == 4) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output POINT_DATA again after intervening CELL_DATA");

339:       if (doOutput) {
340:         PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", n/bs);
341:       }
342:     } else {
343:       if (outputState == 0) {
344:         outputState = 2;
345:         doOutput = 1;
346:       } else if (outputState == 1) {
347:         outputState = 4;
348:         doOutput = 1;
349:       } else if (outputState == 2) doOutput = 0;
350:       else if (outputState == 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Tried to output CELL_DATA again after intervening POINT_DATA");
351:       else if (outputState == 4) doOutput = 0;

353:       if (doOutput) {
354:         PetscViewerASCIIPrintf(viewer, "CELL_DATA %d\n", n);
355:       }
356:     }
357:     PetscObjectComposedDataSetInt((PetscObject) viewer, stateId, outputState);
358:     if (name) {
359:       if (bs == 3) {
360:         PetscViewerASCIIPrintf(viewer, "VECTORS %s double\n", name);
361:       } else {
362:         PetscViewerASCIIPrintf(viewer, "SCALARS %s double %d\n", name, bs);
363:       }
364:     } else {
365:       PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", bs);
366:     }
367:     if (bs != 3) {
368:       PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
369:     }
370:     for (i=0; i<n/bs; i++) {
371:       for (b=0; b<bs; b++) {
372:         if (b > 0) {
373:           PetscViewerASCIIPrintf(viewer," ");
374:         }
375: #if !defined(PETSC_USE_COMPLEX)
376:         PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
377: #endif
378:       }
379:       PetscViewerASCIIPrintf(viewer,"\n");
380:     }
381:   } else if (format == PETSC_VIEWER_ASCII_VTK_COORDS) {
382:     PetscInt bs, b;

384:     VecGetBlockSize(xin, &bs);
385:     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);
386:     for (i=0; i<n/bs; i++) {
387:       for (b=0; b<bs; b++) {
388:         if (b > 0) {
389:           PetscViewerASCIIPrintf(viewer," ");
390:         }
391: #if !defined(PETSC_USE_COMPLEX)
392:         PetscViewerASCIIPrintf(viewer,"%g",(double)xv[i*bs+b]);
393: #endif
394:       }
395:       for (b=bs; b<3; b++) {
396:         PetscViewerASCIIPrintf(viewer," 0.0");
397:       }
398:       PetscViewerASCIIPrintf(viewer,"\n");
399:     }
400:   } else if (format == PETSC_VIEWER_ASCII_PCICE) {
401:     PetscInt bs, b;

403:     VecGetBlockSize(xin, &bs);
404:     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);
405:     PetscViewerASCIIPrintf(viewer,"%D\n", xin->map->N/bs);
406:     for (i=0; i<n/bs; i++) {
407:       PetscViewerASCIIPrintf(viewer,"%7D   ", i+1);
408:       for (b=0; b<bs; b++) {
409:         if (b > 0) {
410:           PetscViewerASCIIPrintf(viewer," ");
411:         }
412: #if !defined(PETSC_USE_COMPLEX)
413:         PetscViewerASCIIPrintf(viewer,"% 12.5E",(double)xv[i*bs+b]);
414: #endif
415:       }
416:       PetscViewerASCIIPrintf(viewer,"\n");
417:     }
418:   } else if (format == PETSC_VIEWER_ASCII_GLVIS) {
419:     /* GLVis ASCII visualization/dump: this function mimicks mfem::GridFunction::Save() */
420:     const PetscScalar       *array;
421:     PetscInt                i,n,vdim, ordering = 1; /* mfem::FiniteElementSpace::Ordering::byVDIM */
422:     PetscContainer          glvis_container;
423:     PetscViewerGLVisVecInfo glvis_vec_info;
424:     PetscViewerGLVisInfo    glvis_info;
425:     PetscErrorCode          ierr;

427:     /* mfem::FiniteElementSpace::Save() */
428:     VecGetBlockSize(xin,&vdim);
429:     PetscViewerASCIIPrintf(viewer,"FiniteElementSpace\n");
430:     PetscObjectQuery((PetscObject)xin,"_glvis_info_container",(PetscObject*)&glvis_container);
431:     if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_PLIB,"Missing GLVis container");
432:     PetscContainerGetPointer(glvis_container,(void**)&glvis_vec_info);
433:     PetscViewerASCIIPrintf(viewer,"%s\n",glvis_vec_info->fec_type);
434:     PetscViewerASCIIPrintf(viewer,"VDim: %d\n",vdim);
435:     PetscViewerASCIIPrintf(viewer,"Ordering: %d\n",ordering);
436:     PetscViewerASCIIPrintf(viewer,"\n");
437:     /* mfem::Vector::Print() */
438:     PetscObjectQuery((PetscObject)viewer,"_glvis_info_container",(PetscObject*)&glvis_container);
439:     if (!glvis_container) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_PLIB,"Missing GLVis container");
440:     PetscContainerGetPointer(glvis_container,(void**)&glvis_info);
441:     if (glvis_info->enabled) {
442:       VecGetLocalSize(xin,&n);
443:       VecGetArrayRead(xin,&array);
444:       for (i=0;i<n;i++) {
445:         PetscViewerASCIIPrintf(viewer,glvis_info->fmt,(double)PetscRealPart(array[i]));
446:         PetscViewerASCIIPrintf(viewer,"\n");
447:       }
448:       VecRestoreArrayRead(xin,&array);
449:     }
450:   } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
451:     /* No info */
452:   } else {
453:     for (i=0; i<n; i++) {
454:       if (format == PETSC_VIEWER_ASCII_INDEX) {
455:         PetscViewerASCIIPrintf(viewer,"%D: ",i);
456:       }
457: #if defined(PETSC_USE_COMPLEX)
458:       if (PetscImaginaryPart(xv[i]) > 0.0) {
459:         PetscViewerASCIIPrintf(viewer,"%g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
460:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
461:         PetscViewerASCIIPrintf(viewer,"%g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
462:       } else {
463:         PetscViewerASCIIPrintf(viewer,"%g\n",(double)PetscRealPart(xv[i]));
464:       }
465: #else
466:       PetscViewerASCIIPrintf(viewer,"%g\n",(double)xv[i]);
467: #endif
468:     }
469:   }
470:   PetscViewerFlush(viewer);
471:   VecRestoreArrayRead(xin,&xv);
472:   return(0);
473: }

475:  #include <petscdraw.h>
476: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
477: {
478:   PetscDraw         draw;
479:   PetscBool         isnull;
480:   PetscDrawLG       lg;
481:   PetscErrorCode    ierr;
482:   PetscInt          i,c,bs = PetscAbs(xin->map->bs),n = xin->map->n/bs;
483:   const PetscScalar *xv;
484:   PetscReal         *xx,*yy,xmin,xmax,h;
485:   int               colors[] = {PETSC_DRAW_RED};
486:   PetscViewerFormat format;
487:   PetscDrawAxis     axis;

490:   PetscViewerDrawGetDraw(v,0,&draw);
491:   PetscDrawIsNull(draw,&isnull);
492:   if (isnull) return(0);

494:   PetscViewerGetFormat(v,&format);
495:   PetscMalloc2(n,&xx,n,&yy);
496:   VecGetArrayRead(xin,&xv);
497:   for (c=0; c<bs; c++) {
498:     PetscViewerDrawGetDrawLG(v,c,&lg);
499:     PetscDrawLGReset(lg);
500:     PetscDrawLGSetDimension(lg,1);
501:     PetscDrawLGSetColors(lg,colors);
502:     if (format == PETSC_VIEWER_DRAW_LG_XRANGE) {
503:       PetscDrawLGGetAxis(lg,&axis);
504:       PetscDrawAxisGetLimits(axis,&xmin,&xmax,NULL,NULL);
505:       h = (xmax - xmin)/n;
506:       for (i=0; i<n; i++) xx[i] = i*h + 0.5*h; /* cell center */
507:     } else {
508:       for (i=0; i<n; i++) xx[i] = (PetscReal)i;
509:     }
510:     for (i=0; i<n; i++) yy[i] = PetscRealPart(xv[c + i*bs]);

512:     PetscDrawLGAddPoints(lg,n,&xx,&yy);
513:     PetscDrawLGDraw(lg);
514:     PetscDrawLGSave(lg);
515:   }
516:   VecRestoreArrayRead(xin,&xv);
517:   PetscFree2(xx,yy);
518:   return(0);
519: }

521: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
522: {
523:   PetscErrorCode    ierr;
524:   PetscDraw         draw;
525:   PetscBool         isnull;

528:   PetscViewerDrawGetDraw(v,0,&draw);
529:   PetscDrawIsNull(draw,&isnull);
530:   if (isnull) return(0);

532:   VecView_Seq_Draw_LG(xin,v);
533:   return(0);
534: }

536: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
537: {
538:   PetscErrorCode    ierr;
539:   int               fdes;
540:   PetscInt          n = xin->map->n,classid=VEC_FILE_CLASSID;
541:   FILE              *file;
542:   const PetscScalar *xv;
543: #if defined(PETSC_HAVE_MPIIO)
544:   PetscBool         isMPIIO;
545: #endif
546:   PetscBool         skipHeader;
547:   PetscViewerFormat format;

550:   /* Write vector header */
551:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
552:   if (!skipHeader) {
553:     PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
554:     PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
555:   }

557:   /* Write vector contents */
558: #if defined(PETSC_HAVE_MPIIO)
559:   PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);
560:   if (!isMPIIO) {
561: #endif
562:     PetscViewerBinaryGetDescriptor(viewer,&fdes);
563:     VecGetArrayRead(xin,&xv);
564:     PetscBinaryWrite(fdes,(void*)xv,n,PETSC_SCALAR,PETSC_FALSE);
565:     VecRestoreArrayRead(xin,&xv);
566:     PetscViewerGetFormat(viewer,&format);
567:     if (format == PETSC_VIEWER_BINARY_MATLAB) {
568:       MPI_Comm   comm;
569:       FILE       *info;
570:       const char *name;

572:       PetscObjectGetName((PetscObject)xin,&name);
573:       PetscObjectGetComm((PetscObject)viewer,&comm);
574:       PetscViewerBinaryGetInfoPointer(viewer,&info);
575:       PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
576:       PetscFPrintf(comm,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
577:       PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
578:     }
579: #if defined(PETSC_HAVE_MPIIO)
580:   } else {
581:     MPI_Offset   off;
582:     MPI_File     mfdes;
583:     PetscMPIInt  lsize;

585:     PetscMPIIntCast(n,&lsize);
586:     PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
587:     PetscViewerBinaryGetMPIIOOffset(viewer,&off);
588:     MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
589:     VecGetArrayRead(xin,&xv);
590:     MPIU_File_write_all(mfdes,(void*)xv,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
591:     VecRestoreArrayRead(xin,&xv);
592:     PetscViewerBinaryAddMPIIOOffset(viewer,n*sizeof(PetscScalar));
593:   }
594: #endif

596:   PetscViewerBinaryGetInfoPointer(viewer,&file);
597:   if (file) {
598:     if (((PetscObject)xin)->prefix) {
599:       PetscFPrintf(PETSC_COMM_SELF,file,"-%svecload_block_size %D\n",((PetscObject)xin)->prefix,PetscAbs(xin->map->bs));
600:     } else {
601:       PetscFPrintf(PETSC_COMM_SELF,file,"-vecload_block_size %D\n",PetscAbs(xin->map->bs));
602:     }
603:   }
604:   return(0);
605: }

607: #if defined(PETSC_HAVE_MATLAB_ENGINE)
608:  #include <petscmatlab.h>
609: #include <mat.h>   /* MATLAB include file */
610: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
611: {
612:   PetscErrorCode    ierr;
613:   PetscInt          n;
614:   const PetscScalar *array;

617:   VecGetLocalSize(vec,&n);
618:   PetscObjectName((PetscObject)vec);
619:   VecGetArrayRead(vec,&array);
620:   PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
621:   VecRestoreArrayRead(vec,&array);
622:   return(0);
623: }
624: #endif

626: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
627: {
629:   PetscBool      isdraw,iascii,issocket,isbinary;
630: #if defined(PETSC_HAVE_MATHEMATICA)
631:   PetscBool      ismathematica;
632: #endif
633: #if defined(PETSC_HAVE_MATLAB_ENGINE)
634:   PetscBool      ismatlab;
635: #endif
636: #if defined(PETSC_HAVE_HDF5)
637:   PetscBool      ishdf5;
638: #endif
639:   PetscBool      isglvis;
640: #if defined(PETSC_HAVE_ADIOS)
641:   PetscBool      isadios;
642: #endif
643: #if defined(PETSC_HAVE_ADIOS2)
644:   PetscBool      isadios2;
645: #endif

648:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
649:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
650:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
651:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
652: #if defined(PETSC_HAVE_MATHEMATICA)
653:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
654: #endif
655: #if defined(PETSC_HAVE_HDF5)
656:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
657: #endif
658: #if defined(PETSC_HAVE_MATLAB_ENGINE)
659:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
660: #endif
661:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
662: #if defined(PETSC_HAVE_ADIOS)
663:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
664: #endif
665: #if defined(PETSC_HAVE_ADIOS2)
666:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS2,&isadios2);
667: #endif

669:   if (isdraw) {
670:     VecView_Seq_Draw(xin,viewer);
671:   } else if (iascii) {
672:     VecView_Seq_ASCII(xin,viewer);
673:   } else if (isbinary) {
674:     VecView_Seq_Binary(xin,viewer);
675: #if defined(PETSC_HAVE_MATHEMATICA)
676:   } else if (ismathematica) {
677:     PetscViewerMathematicaPutVector(viewer,xin);
678: #endif
679: #if defined(PETSC_HAVE_HDF5)
680:   } else if (ishdf5) {
681:     VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
682: #endif
683: #if defined(PETSC_HAVE_ADIOS)
684:   } else if (isadios) {
685:     VecView_MPI_ADIOS(xin,viewer); /* Reusing VecView_MPI_ADIOS ... don't want code duplication*/
686: #endif
687: #if defined(PETSC_HAVE_ADIOS2)
688:   } else if (isadios2) {
689:     VecView_MPI_ADIOS2(xin,viewer); /* Reusing VecView_MPI_ADIOS2 ... don't want code duplication*/
690: #endif
691: #if defined(PETSC_HAVE_MATLAB_ENGINE)
692:   } else if (ismatlab) {
693:     VecView_Seq_Matlab(xin,viewer);
694: #endif
695:   } else if (isglvis) {
696:     VecView_GLVis(xin,viewer);
697:   }
698:   return(0);
699: }

701: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
702: {
703:   const PetscScalar *xx;
704:   PetscInt          i;
705:   PetscErrorCode    ierr;

708:   VecGetArrayRead(xin,&xx);
709:   for (i=0; i<ni; i++) {
710:     if (xin->stash.ignorenegidx && ix[i] < 0) continue;
711: #if defined(PETSC_USE_DEBUG)
712:     if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
713:     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);
714: #endif
715:     y[i] = xx[ix[i]];
716:   }
717:   VecRestoreArrayRead(xin,&xx);
718:   return(0);
719: }

721: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
722: {
723:   PetscScalar    *xx;
724:   PetscInt       i;

728:   VecGetArray(xin,&xx);
729:   if (m == INSERT_VALUES) {
730:     for (i=0; i<ni; i++) {
731:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
732: #if defined(PETSC_USE_DEBUG)
733:       if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
734:       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);
735: #endif
736:       xx[ix[i]] = y[i];
737:     }
738:   } else {
739:     for (i=0; i<ni; i++) {
740:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
741: #if defined(PETSC_USE_DEBUG)
742:       if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
743:       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);
744: #endif
745:       xx[ix[i]] += y[i];
746:     }
747:   }
748:   VecRestoreArray(xin,&xx);
749:   return(0);
750: }

752: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
753: {
754:   PetscScalar    *xx,*y = (PetscScalar*)yin;
755:   PetscInt       i,bs,start,j;

758:   /*
759:        For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
760:   */
762:   VecGetBlockSize(xin,&bs);
763:   VecGetArray(xin,&xx);
764:   if (m == INSERT_VALUES) {
765:     for (i=0; i<ni; i++, y+=bs) {
766:       start = bs*ix[i];
767:       if (start < 0) continue;
768: #if defined(PETSC_USE_DEBUG)
769:       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);
770: #endif
771:       for (j=0; j<bs; j++) xx[start+j] = y[j];
772:     }
773:   } else {
774:     for (i=0; i<ni; i++, y+=bs) {
775:       start = bs*ix[i];
776:       if (start < 0) continue;
777: #if defined(PETSC_USE_DEBUG)
778:       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);
779: #endif
780:       for (j=0; j<bs; j++) xx[start+j] += y[j];
781:     }
782:   }
783:   VecRestoreArray(xin,&xx);
784:   return(0);
785: }

787: PetscErrorCode VecDestroy_Seq(Vec v)
788: {
789:   Vec_Seq        *vs = (Vec_Seq*)v->data;

793: #if defined(PETSC_USE_LOG)
794:   PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
795: #endif
796:   PetscFree(vs->array_allocated);
797:   PetscFree(v->data);
798:   return(0);
799: }

801: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
802: {
804:   if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
805:   return(0);
806: }

808: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
809: {

813:   VecCreate(PetscObjectComm((PetscObject)win),V);
814:   VecSetSizes(*V,win->map->n,win->map->n);
815:   VecSetType(*V,((PetscObject)win)->type_name);
816:   PetscLayoutReference(win->map,&(*V)->map);
817:   PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
818:   PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);

820:   (*V)->ops->view          = win->ops->view;
821:   (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
822:   return(0);
823: }

825: static struct _VecOps DvOps = {VecDuplicate_Seq, /* 1 */
826:                                VecDuplicateVecs_Default,
827:                                VecDestroyVecs_Default,
828:                                VecDot_Seq,
829:                                VecMDot_Seq,
830:                                VecNorm_Seq,
831:                                VecTDot_Seq,
832:                                VecMTDot_Seq,
833:                                VecScale_Seq,
834:                                VecCopy_Seq, /* 10 */
835:                                VecSet_Seq,
836:                                VecSwap_Seq,
837:                                VecAXPY_Seq,
838:                                VecAXPBY_Seq,
839:                                VecMAXPY_Seq,
840:                                VecAYPX_Seq,
841:                                VecWAXPY_Seq,
842:                                VecAXPBYPCZ_Seq,
843:                                VecPointwiseMult_Seq,
844:                                VecPointwiseDivide_Seq,
845:                                VecSetValues_Seq, /* 20 */
846:                                0,0,
847:                                0,
848:                                VecGetSize_Seq,
849:                                VecGetSize_Seq,
850:                                0,
851:                                VecMax_Seq,
852:                                VecMin_Seq,
853:                                VecSetRandom_Seq,
854:                                VecSetOption_Seq, /* 30 */
855:                                VecSetValuesBlocked_Seq,
856:                                VecDestroy_Seq,
857:                                VecView_Seq,
858:                                VecPlaceArray_Seq,
859:                                VecReplaceArray_Seq,
860:                                VecDot_Seq,
861:                                VecTDot_Seq,
862:                                VecNorm_Seq,
863:                                VecMDot_Seq,
864:                                VecMTDot_Seq, /* 40 */
865:                                VecLoad_Default,
866:                                VecReciprocal_Default,
867:                                VecConjugate_Seq,
868:                                0,
869:                                0,
870:                                VecResetArray_Seq,
871:                                0,
872:                                VecMaxPointwiseDivide_Seq,
873:                                VecPointwiseMax_Seq,
874:                                VecPointwiseMaxAbs_Seq,
875:                                VecPointwiseMin_Seq,
876:                                VecGetValues_Seq,
877:                                0,
878:                                0,
879:                                0,
880:                                0,
881:                                0,
882:                                0,
883:                                VecStrideGather_Default,
884:                                VecStrideScatter_Default,
885:                                0,
886:                                0,
887:                                0,
888:                                0,
889:                                0,
890:                                VecStrideSubSetGather_Default,
891:                                VecStrideSubSetScatter_Default,
892:                                0,
893:                                0
894: };


897: /*
898:       This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
899: */
900: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
901: {
902:   Vec_Seq        *s;

906:   PetscNewLog(v,&s);
907:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));

909:   v->data            = (void*)s;
910:   v->petscnative     = PETSC_TRUE;
911:   s->array           = (PetscScalar*)array;
912:   s->array_allocated = 0;

914:   PetscLayoutSetUp(v->map);
915:   PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
916: #if defined(PETSC_HAVE_MATLAB_ENGINE)
917:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
918:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
919: #endif
920:   return(0);
921: }

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

927:    Collective

929:    Input Parameter:
930: +  comm - the communicator, should be PETSC_COMM_SELF
931: .  bs - the block size
932: .  n - the vector length
933: -  array - memory where the vector elements are to be stored.

935:    Output Parameter:
936: .  V - the vector

938:    Notes:
939:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
940:    same type as an existing vector.

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

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

948:    Level: intermediate

950: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
951:           VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
952: @*/
953: PetscErrorCode  VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
954: {
956:   PetscMPIInt    size;

959:   VecCreate(comm,V);
960:   VecSetSizes(*V,n,n);
961:   VecSetBlockSize(*V,bs);
962:   MPI_Comm_size(comm,&size);
963:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
964:   VecCreate_Seq_Private(*V,array);
965:   return(0);
966: }