Actual source code: bvec2.c

petsc-3.8.3 2017-12-09
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++) ww[i] = xx[i] / yy[i];

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

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

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

139: PetscErrorCode VecGetSize_Seq(Vec vin,PetscInt *size)
140: {
142:   *size = vin->map->n;
143:   return(0);
144: }



148: PetscErrorCode VecConjugate_Seq(Vec xin)
149: {
150:   PetscScalar    *x;
151:   PetscInt       n = xin->map->n;

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

164: PetscErrorCode VecResetArray_Seq(Vec vin)
165: {
166:   Vec_Seq *v = (Vec_Seq*)vin->data;

169:   v->array         = v->unplacedarray;
170:   v->unplacedarray = 0;
171:   return(0);
172: }

174: PetscErrorCode VecCopy_Seq(Vec xin,Vec yin)
175: {
176:   PetscScalar       *ya;
177:   const PetscScalar *xa;
178:   PetscErrorCode    ierr;

181:   if (xin != yin) {
182:     VecGetArrayRead(xin,&xa);
183:     VecGetArray(yin,&ya);
184:     PetscMemcpy(ya,xa,xin->map->n*sizeof(PetscScalar));
185:     VecRestoreArrayRead(xin,&xa);
186:     VecRestoreArray(yin,&ya);
187:   }
188:   return(0);
189: }

191: PetscErrorCode VecSwap_Seq(Vec xin,Vec yin)
192: {
193:   PetscScalar    *ya, *xa;
195:   PetscBLASInt   one = 1,bn;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

485:   PetscViewerDrawGetDraw(v,0,&draw);
486:   PetscDrawIsNull(draw,&isnull);
487:   if (isnull) return(0);

489:   PetscMalloc2(n,&xx,n,&yy);
490:   VecGetArrayRead(xin,&xv);
491:   for (c=0; c<bs; c++) {
492:     PetscViewerDrawGetDrawLG(v,c,&lg);
493:     PetscDrawLGReset(lg);
494:     PetscDrawLGSetDimension(lg,1);
495:     PetscDrawLGSetColors(lg,colors);
496:     for (i=0; i<n; i++) {
497:       xx[i] = (PetscReal)i;
498:       yy[i] = PetscRealPart(xv[c + i*bs]);
499:     }
500:     PetscDrawLGAddPoints(lg,n,&xx,&yy);
501:     PetscDrawLGDraw(lg);
502:     PetscDrawLGSave(lg);
503:   }
504:   VecRestoreArrayRead(xin,&xv);
505:   PetscFree2(xx,yy);
506:   return(0);
507: }

509: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
510: {
511:   PetscErrorCode    ierr;
512:   PetscDraw         draw;
513:   PetscBool         isnull;

516:   PetscViewerDrawGetDraw(v,0,&draw);
517:   PetscDrawIsNull(draw,&isnull);
518:   if (isnull) return(0);
519:   PetscViewerPushFormat(v,PETSC_VIEWER_DRAW_LG);
520:   VecView_Seq_Draw_LG(xin,v);
521:   PetscViewerPopFormat(v);
522:   return(0);
523: }

525: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
526: {
527:   PetscErrorCode    ierr;
528:   int               fdes;
529:   PetscInt          n = xin->map->n,classid=VEC_FILE_CLASSID;
530:   FILE              *file;
531:   const PetscScalar *xv;
532: #if defined(PETSC_HAVE_MPIIO)
533:   PetscBool         isMPIIO;
534: #endif
535:   PetscBool         skipHeader;
536:   PetscViewerFormat format;

539:   /* Write vector header */
540:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
541:   if (!skipHeader) {
542:     PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
543:     PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
544:   }

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

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

574:     PetscMPIIntCast(n,&lsize);
575:     PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
576:     PetscViewerBinaryGetMPIIOOffset(viewer,&off);
577:     MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
578:     VecGetArrayRead(xin,&xv);
579:     MPIU_File_write_all(mfdes,(void*)xv,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
580:     VecRestoreArrayRead(xin,&xv);
581:     PetscViewerBinaryAddMPIIOOffset(viewer,n*sizeof(PetscScalar));
582:   }
583: #endif

585:   PetscViewerBinaryGetInfoPointer(viewer,&file);
586:   if (file) {
587:     if (((PetscObject)xin)->prefix) {
588:       PetscFPrintf(PETSC_COMM_SELF,file,"-%svecload_block_size %D\n",((PetscObject)xin)->prefix,PetscAbs(xin->map->bs));
589:     } else {
590:       PetscFPrintf(PETSC_COMM_SELF,file,"-vecload_block_size %D\n",PetscAbs(xin->map->bs));
591:     }
592:   }
593:   return(0);
594: }

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

606:   VecGetLocalSize(vec,&n);
607:   PetscObjectName((PetscObject)vec);
608:   VecGetArrayRead(vec,&array);
609:   PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
610:   VecRestoreArrayRead(vec,&array);
611:   return(0);
612: }
613: #endif

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

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

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

670: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
671: {
672:   const PetscScalar *xx;
673:   PetscInt          i;
674:   PetscErrorCode    ierr;

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

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

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

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: }

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

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

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

779: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
780: {

784:   VecCreate(PetscObjectComm((PetscObject)win),V);
785:   VecSetSizes(*V,win->map->n,win->map->n);
786:   VecSetType(*V,((PetscObject)win)->type_name);
787:   PetscLayoutReference(win->map,&(*V)->map);
788:   PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
789:   PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);

791:   (*V)->ops->view          = win->ops->view;
792:   (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
793:   return(0);
794: }

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


868: /*
869:       This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
870: */
871: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
872: {
873:   Vec_Seq        *s;

877:   PetscNewLog(v,&s);
878:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));

880:   v->data            = (void*)s;
881:   v->petscnative     = PETSC_TRUE;
882:   s->array           = (PetscScalar*)array;
883:   s->array_allocated = 0;

885:   PetscLayoutSetUp(v->map);
886:   PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
887: #if defined(PETSC_HAVE_MATLAB_ENGINE)
888:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
889:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
890: #endif
891:   return(0);
892: }

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

898:    Collective on MPI_Comm

900:    Input Parameter:
901: +  comm - the communicator, should be PETSC_COMM_SELF
902: .  bs - the block size
903: .  n - the vector length
904: -  array - memory where the vector elements are to be stored.

906:    Output Parameter:
907: .  V - the vector

909:    Notes:
910:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
911:    same type as an existing vector.

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

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

919:    Level: intermediate

921:    Concepts: vectors^creating with array

923: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
924:           VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
925: @*/
926: PetscErrorCode  VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
927: {
929:   PetscMPIInt    size;

932:   VecCreate(comm,V);
933:   VecSetSizes(*V,n,n);
934:   VecSetBlockSize(*V,bs);
935:   MPI_Comm_size(comm,&size);
936:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
937:   VecCreate_Seq_Private(*V,array);
938:   return(0);
939: }