Actual source code: bvec2.c

petsc-3.6.1 2015-07-22
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:     VecGetArrayRead(xin,&xx);
263:     PetscStackCallBLAS("BLASasum",*z   = BLASasum_(&bn,xx,&one));
264:     VecRestoreArrayRead(xin,&xx);
265:     PetscLogFlops(PetscMax(n-1.0,0.0));
266:   } else if (type == NORM_1_AND_2) {
267:     VecNorm_Seq(xin,NORM_1,z);
268:     VecNorm_Seq(xin,NORM_2,z+1);
269:   }
270:   return(0);
271: }

275: PetscErrorCode VecView_Seq_ASCII(Vec xin,PetscViewer viewer)
276: {
277:   PetscErrorCode    ierr;
278:   PetscInt          i,n = xin->map->n;
279:   const char        *name;
280:   PetscViewerFormat format;
281:   const PetscScalar *xv;

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

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

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

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

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

408:     VecGetBlockSize(xin, &bs);
409:     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);
410:     PetscViewerASCIIPrintf(viewer,"%D\n", xin->map->N/bs);
411:     for (i=0; i<n/bs; i++) {
412:       PetscViewerASCIIPrintf(viewer,"%7D   ", i+1);
413:       for (b=0; b<bs; b++) {
414:         if (b > 0) {
415:           PetscViewerASCIIPrintf(viewer," ");
416:         }
417: #if !defined(PETSC_USE_COMPLEX)
418:         PetscViewerASCIIPrintf(viewer,"% 12.5E",(double)xv[i*bs+b]);
419: #endif
420:       }
421:       PetscViewerASCIIPrintf(viewer,"\n");
422:     }
423:   } else if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
424:     /* No info */
425:   } else {
426:     for (i=0; i<n; i++) {
427:       if (format == PETSC_VIEWER_ASCII_INDEX) {
428:         PetscViewerASCIIPrintf(viewer,"%D: ",i);
429:       }
430: #if defined(PETSC_USE_COMPLEX)
431:       if (PetscImaginaryPart(xv[i]) > 0.0) {
432:         PetscViewerASCIIPrintf(viewer,"%g + %g i\n",(double)PetscRealPart(xv[i]),(double)PetscImaginaryPart(xv[i]));
433:       } else if (PetscImaginaryPart(xv[i]) < 0.0) {
434:         PetscViewerASCIIPrintf(viewer,"%g - %g i\n",(double)PetscRealPart(xv[i]),-(double)PetscImaginaryPart(xv[i]));
435:       } else {
436:         PetscViewerASCIIPrintf(viewer,"%g\n",(double)PetscRealPart(xv[i]));
437:       }
438: #else
439:       PetscViewerASCIIPrintf(viewer,"%g\n",(double)xv[i]);
440: #endif
441:     }
442:   }
443:   PetscViewerFlush(viewer);
444:   VecRestoreArrayRead(xin,&xv);
445:   return(0);
446: }

448: #include <petscdraw.h>
451: PetscErrorCode VecView_Seq_Draw_LG(Vec xin,PetscViewer v)
452: {
453:   PetscErrorCode    ierr;
454:   PetscInt          i,c,bs = PetscAbs(xin->map->bs),n = xin->map->n/bs;
455:   PetscDraw         win;
456:   PetscReal         *xx;
457:   PetscDrawLG       lg;
458:   const PetscScalar *xv;
459:   PetscReal         *yy;

462:   PetscMalloc1(n,&xx);
463:   PetscMalloc1(n,&yy);
464:   VecGetArrayRead(xin,&xv);
465:   for (c=0; c<bs; c++) {
466:     PetscViewerDrawGetDrawLG(v,c,&lg);
467:     PetscDrawLGGetDraw(lg,&win);
468:     PetscDrawCheckResizedWindow(win);
469:     PetscDrawLGReset(lg);
470:     for (i=0; i<n; i++) {
471:       xx[i] = (PetscReal) i;
472:       yy[i] = PetscRealPart(xv[c + i*bs]);
473:     }
474:     PetscDrawLGAddPoints(lg,n,&xx,&yy);
475:     PetscDrawLGDraw(lg);
476:     PetscDrawSynchronizedFlush(win);
477:   }
478:   VecRestoreArrayRead(xin,&xv);
479:   PetscFree(yy);
480:   PetscFree(xx);
481:   return(0);
482: }

486: PetscErrorCode VecView_Seq_Draw(Vec xin,PetscViewer v)
487: {
488:   PetscErrorCode    ierr;
489:   PetscDraw         draw;
490:   PetscBool         isnull;
491:   PetscViewerFormat format;

494:   PetscViewerDrawGetDraw(v,0,&draw);
495:   PetscDrawIsNull(draw,&isnull); if (isnull) return(0);

497:   PetscViewerGetFormat(v,&format);
498:   /*
499:      Currently it only supports drawing to a line graph */
500:   if (format != PETSC_VIEWER_DRAW_LG) {
501:     PetscViewerPushFormat(v,PETSC_VIEWER_DRAW_LG);
502:   }
503:   VecView_Seq_Draw_LG(xin,v);
504:   if (format != PETSC_VIEWER_DRAW_LG) {
505:     PetscViewerPopFormat(v);
506:   }
507:   return(0);
508: }

512: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
513: {
514:   PetscErrorCode    ierr;
515:   int               fdes;
516:   PetscInt          n = xin->map->n,classid=VEC_FILE_CLASSID;
517:   FILE              *file;
518:   const PetscScalar *xv;
519: #if defined(PETSC_HAVE_MPIIO)
520:   PetscBool         isMPIIO;
521: #endif
522:   PetscBool         skipHeader;
523:   PetscViewerFormat format;

526:   /* Write vector header */
527:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
528:   if (!skipHeader) {
529:     PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
530:     PetscViewerBinaryWrite(viewer,&n,1,PETSC_INT,PETSC_FALSE);
531:   }

533:   /* Write vector contents */
534: #if defined(PETSC_HAVE_MPIIO)
535:   PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);
536:   if (!isMPIIO) {
537: #endif
538:     PetscViewerBinaryGetDescriptor(viewer,&fdes);
539:     VecGetArrayRead(xin,&xv);
540:     PetscBinaryWrite(fdes,(void*)xv,n,PETSC_SCALAR,PETSC_FALSE);
541:     VecRestoreArrayRead(xin,&xv);
542:     PetscViewerGetFormat(viewer,&format);
543:     if (format == PETSC_VIEWER_BINARY_MATLAB) {
544:       MPI_Comm   comm;
545:       FILE       *info;
546:       const char *name;

548:       PetscObjectGetName((PetscObject)xin,&name);
549:       PetscObjectGetComm((PetscObject)viewer,&comm);
550:       PetscViewerBinaryGetInfoPointer(viewer,&info);
551:       PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
552:       PetscFPrintf(comm,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
553:       PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
554:     }
555: #if defined(PETSC_HAVE_MPIIO)
556:   } else {
557:     MPI_Offset   off;
558:     MPI_File     mfdes;
559:     PetscMPIInt  lsize;

561:     PetscMPIIntCast(n,&lsize);
562:     PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
563:     PetscViewerBinaryGetMPIIOOffset(viewer,&off);
564:     MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
565:     VecGetArrayRead(xin,&xv);
566:     MPIU_File_write_all(mfdes,(void*)xv,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
567:     VecRestoreArrayRead(xin,&xv);
568:     PetscViewerBinaryAddMPIIOOffset(viewer,n*sizeof(PetscScalar));
569:   }
570: #endif

572:   PetscViewerBinaryGetInfoPointer(viewer,&file);
573:   if (file) {
574:     if (((PetscObject)xin)->prefix) {
575:       PetscFPrintf(PETSC_COMM_SELF,file,"-%svecload_block_size %D\n",((PetscObject)xin)->prefix,PetscAbs(xin->map->bs));
576:     } else {
577:       PetscFPrintf(PETSC_COMM_SELF,file,"-vecload_block_size %D\n",PetscAbs(xin->map->bs));
578:     }
579:   }
580:   return(0);
581: }

583: #if defined(PETSC_HAVE_MATLAB_ENGINE)
584: #include <petscmatlab.h>
585: #include <mat.h>   /* MATLAB include file */
588: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
589: {
590:   PetscErrorCode    ierr;
591:   PetscInt          n;
592:   const PetscScalar *array;

595:   VecGetLocalSize(vec,&n);
596:   PetscObjectName((PetscObject)vec);
597:   VecGetArrayRead(vec,&array);
598:   PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
599:   VecRestoreArrayRead(vec,&array);
600:   return(0);
601: }
602: #endif

606: PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
607: {
609:   PetscBool      isdraw,iascii,issocket,isbinary;
610: #if defined(PETSC_HAVE_MATHEMATICA)
611:   PetscBool      ismathematica;
612: #endif
613: #if defined(PETSC_HAVE_MATLAB_ENGINE)
614:   PetscBool      ismatlab;
615: #endif
616: #if defined(PETSC_HAVE_HDF5)
617:   PetscBool      ishdf5;
618: #endif

621:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
622:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
623:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
624:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
625: #if defined(PETSC_HAVE_MATHEMATICA)
626:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
627: #endif
628: #if defined(PETSC_HAVE_HDF5)
629:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
630: #endif
631: #if defined(PETSC_HAVE_MATLAB_ENGINE)
632:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
633: #endif

635:   if (isdraw) {
636:     VecView_Seq_Draw(xin,viewer);
637:   } else if (iascii) {
638:     VecView_Seq_ASCII(xin,viewer);
639:   } else if (isbinary) {
640:     VecView_Seq_Binary(xin,viewer);
641: #if defined(PETSC_HAVE_MATHEMATICA)
642:   } else if (ismathematica) {
643:     PetscViewerMathematicaPutVector(viewer,xin);
644: #endif
645: #if defined(PETSC_HAVE_HDF5)
646:   } else if (ishdf5) {
647:     VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
648: #endif
649: #if defined(PETSC_HAVE_MATLAB_ENGINE)
650:   } else if (ismatlab) {
651:     VecView_Seq_Matlab(xin,viewer);
652: #endif
653:   }
654:   return(0);
655: }

659: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
660: {
661:   const PetscScalar *xx;
662:   PetscInt          i;
663:   PetscErrorCode    ierr;

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

681: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
682: {
683:   PetscScalar    *xx;
684:   PetscInt       i;

688:   VecGetArray(xin,&xx);
689:   if (m == INSERT_VALUES) {
690:     for (i=0; i<ni; i++) {
691:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
692: #if defined(PETSC_USE_DEBUG)
693:       if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
694:       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);
695: #endif
696:       xx[ix[i]] = y[i];
697:     }
698:   } else {
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:   }
708:   VecRestoreArray(xin,&xx);
709:   return(0);
710: }

714: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
715: {
716:   PetscScalar    *xx,*y = (PetscScalar*)yin;
717:   PetscInt       i,bs,start,j;

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

753: PetscErrorCode VecDestroy_Seq(Vec v)
754: {
755:   Vec_Seq        *vs = (Vec_Seq*)v->data;

759: #if defined(PETSC_USE_LOG)
760:   PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
761: #endif
762:   PetscFree(vs->array_allocated);
763:   PetscFree(v->data);
764:   return(0);
765: }

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

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

783:   VecCreate(PetscObjectComm((PetscObject)win),V);
784:   PetscObjectSetPrecision((PetscObject)*V,((PetscObject)win)->precision);
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: */
873: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
874: {
875:   Vec_Seq        *s;

879:   PetscNewLog(v,&s);
880:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));

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

887:   PetscLayoutSetUp(v->map);
888:   PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
889: #if defined(PETSC_HAVE_MATLAB_ENGINE)
890:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
891:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
892: #endif
893: #if defined(PETSC_USE_MIXED_PRECISION)
894:   ((PetscObject)v)->precision = (PetscPrecision)sizeof(PetscReal);
895: #endif
896:   return(0);
897: }

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

905:    Collective on MPI_Comm

907:    Input Parameter:
908: +  comm - the communicator, should be PETSC_COMM_SELF
909: .  bs - the block size
910: .  n - the vector length
911: -  array - memory where the vector elements are to be stored.

913:    Output Parameter:
914: .  V - the vector

916:    Notes:
917:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
918:    same type as an existing vector.

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

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

926:    Level: intermediate

928:    Concepts: vectors^creating with array

930: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
931:           VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
932: @*/
933: PetscErrorCode  VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
934: {
936:   PetscMPIInt    size;

939:   VecCreate(comm,V);
940:   VecSetSizes(*V,n,n);
941:   VecSetBlockSize(*V,bs);
942:   MPI_Comm_size(comm,&size);
943:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
944:   VecCreate_Seq_Private(*V,array);
945:   return(0);
946: }