Actual source code: bvec2.c

petsc-main 2021-04-20
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>
  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:   VecGetArrayWrite(xin,&xx);
137:   for (i=0; i<n; i++) {PetscRandomGetValue(r,&xx[i]);}
138:   VecRestoreArrayWrite(xin,&xx);
139:   return(0);
140: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

488:   PetscViewerDrawGetDraw(v,0,&draw);
489:   PetscDrawIsNull(draw,&isnull);
490:   if (isnull) return(0);

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

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

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

526:   PetscViewerDrawGetDraw(v,0,&draw);
527:   PetscDrawIsNull(draw,&isnull);
528:   if (isnull) return(0);

530:   VecView_Seq_Draw_LG(xin,v);
531:   return(0);
532: }

534: PetscErrorCode VecView_Seq_Binary(Vec xin,PetscViewer viewer)
535: {
536:   return VecView_Binary(xin,viewer);
537: }

539: #if defined(PETSC_HAVE_MATLAB_ENGINE)
540: #include <petscmatlab.h>
541: #include <mat.h>   /* MATLAB include file */
542: PetscErrorCode VecView_Seq_Matlab(Vec vec,PetscViewer viewer)
543: {
544:   PetscErrorCode    ierr;
545:   PetscInt          n;
546:   const PetscScalar *array;

549:   VecGetLocalSize(vec,&n);
550:   PetscObjectName((PetscObject)vec);
551:   VecGetArrayRead(vec,&array);
552:   PetscViewerMatlabPutArray(viewer,n,1,array,((PetscObject)vec)->name);
553:   VecRestoreArrayRead(vec,&array);
554:   return(0);
555: }
556: #endif

558: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec xin,PetscViewer viewer)
559: {
561:   PetscBool      isdraw,iascii,issocket,isbinary;
562: #if defined(PETSC_HAVE_MATHEMATICA)
563:   PetscBool      ismathematica;
564: #endif
565: #if defined(PETSC_HAVE_MATLAB_ENGINE)
566:   PetscBool      ismatlab;
567: #endif
568: #if defined(PETSC_HAVE_HDF5)
569:   PetscBool      ishdf5;
570: #endif
571:   PetscBool      isglvis;
572: #if defined(PETSC_HAVE_ADIOS)
573:   PetscBool      isadios;
574: #endif

577:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
578:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
579:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);
580:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
581: #if defined(PETSC_HAVE_MATHEMATICA)
582:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATHEMATICA,&ismathematica);
583: #endif
584: #if defined(PETSC_HAVE_HDF5)
585:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
586: #endif
587: #if defined(PETSC_HAVE_MATLAB_ENGINE)
588:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERMATLAB,&ismatlab);
589: #endif
590:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
591: #if defined(PETSC_HAVE_ADIOS)
592:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
593: #endif

595:   if (isdraw) {
596:     VecView_Seq_Draw(xin,viewer);
597:   } else if (iascii) {
598:     VecView_Seq_ASCII(xin,viewer);
599:   } else if (isbinary) {
600:     VecView_Seq_Binary(xin,viewer);
601: #if defined(PETSC_HAVE_MATHEMATICA)
602:   } else if (ismathematica) {
603:     PetscViewerMathematicaPutVector(viewer,xin);
604: #endif
605: #if defined(PETSC_HAVE_HDF5)
606:   } else if (ishdf5) {
607:     VecView_MPI_HDF5(xin,viewer); /* Reusing VecView_MPI_HDF5 ... don't want code duplication*/
608: #endif
609: #if defined(PETSC_HAVE_ADIOS)
610:   } else if (isadios) {
611:     VecView_MPI_ADIOS(xin,viewer); /* Reusing VecView_MPI_ADIOS ... don't want code duplication*/
612: #endif
613: #if defined(PETSC_HAVE_MATLAB_ENGINE)
614:   } else if (ismatlab) {
615:     VecView_Seq_Matlab(xin,viewer);
616: #endif
617:   } else if (isglvis) {
618:     VecView_GLVis(xin,viewer);
619:   }
620:   return(0);
621: }

623: PetscErrorCode VecGetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],PetscScalar y[])
624: {
625:   const PetscScalar *xx;
626:   PetscInt          i;
627:   PetscErrorCode    ierr;

630:   VecGetArrayRead(xin,&xx);
631:   for (i=0; i<ni; i++) {
632:     if (xin->stash.ignorenegidx && ix[i] < 0) continue;
633:     if (PetscDefined(USE_DEBUG)) {
634:       if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
635:       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);
636:     }
637:     y[i] = xx[ix[i]];
638:   }
639:   VecRestoreArrayRead(xin,&xx);
640:   return(0);
641: }

643: PetscErrorCode VecSetValues_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode m)
644: {
645:   PetscScalar    *xx;
646:   PetscInt       i;

650:   VecGetArray(xin,&xx);
651:   if (m == INSERT_VALUES) {
652:     for (i=0; i<ni; i++) {
653:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
654:       if (PetscDefined(USE_DEBUG)) {
655:         if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
656:         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);
657:       }
658:       xx[ix[i]] = y[i];
659:     }
660:   } else {
661:     for (i=0; i<ni; i++) {
662:       if (xin->stash.ignorenegidx && ix[i] < 0) continue;
663:       if (PetscDefined(USE_DEBUG)) {
664:         if (ix[i] < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D cannot be negative",ix[i]);
665:         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);
666:       }
667:       xx[ix[i]] += y[i];
668:     }
669:   }
670:   VecRestoreArray(xin,&xx);
671:   return(0);
672: }

674: PetscErrorCode VecSetValuesBlocked_Seq(Vec xin,PetscInt ni,const PetscInt ix[],const PetscScalar yin[],InsertMode m)
675: {
676:   PetscScalar    *xx,*y = (PetscScalar*)yin;
677:   PetscInt       i,bs,start,j;

680:   /*
681:        For optimization could treat bs = 2, 3, 4, 5 as special cases with loop unrolling
682:   */
684:   VecGetBlockSize(xin,&bs);
685:   VecGetArray(xin,&xx);
686:   if (m == INSERT_VALUES) {
687:     for (i=0; i<ni; i++, y+=bs) {
688:       start = bs*ix[i];
689:       if (start < 0) continue;
690:       if (PetscUnlikelyDebug(start >= xin->map->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",start,xin->map->n);
691:       for (j=0; j<bs; j++) xx[start+j] = y[j];
692:     }
693:   } else {
694:     for (i=0; i<ni; i++, y+=bs) {
695:       start = bs*ix[i];
696:       if (start < 0) continue;
697:       if (PetscUnlikelyDebug(start >= xin->map->n)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Out of range index value %D maximum %D",start,xin->map->n);
698:       for (j=0; j<bs; j++) xx[start+j] += y[j];
699:     }
700:   }
701:   VecRestoreArray(xin,&xx);
702:   return(0);
703: }

705: PetscErrorCode VecDestroy_Seq(Vec v)
706: {
707:   Vec_Seq        *vs = (Vec_Seq*)v->data;

711: #if defined(PETSC_USE_LOG)
712:   PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
713: #endif
714:   if (vs) { PetscFree(vs->array_allocated); }
715:   PetscFree(v->data);
716:   return(0);
717: }

719: PetscErrorCode VecSetOption_Seq(Vec v,VecOption op,PetscBool flag)
720: {
722:   if (op == VEC_IGNORE_NEGATIVE_INDICES) v->stash.ignorenegidx = flag;
723:   return(0);
724: }

726: PetscErrorCode VecDuplicate_Seq(Vec win,Vec *V)
727: {

731:   VecCreate(PetscObjectComm((PetscObject)win),V);
732:   VecSetSizes(*V,win->map->n,win->map->n);
733:   VecSetType(*V,((PetscObject)win)->type_name);
734:   PetscLayoutReference(win->map,&(*V)->map);
735:   PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*V))->olist);
736:   PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*V))->qlist);

738:   (*V)->ops->view          = win->ops->view;
739:   (*V)->stash.ignorenegidx = win->stash.ignorenegidx;
740:   return(0);
741: }

743: static struct _VecOps DvOps = {VecDuplicate_Seq, /* 1 */
744:                                VecDuplicateVecs_Default,
745:                                VecDestroyVecs_Default,
746:                                VecDot_Seq,
747:                                VecMDot_Seq,
748:                                VecNorm_Seq,
749:                                VecTDot_Seq,
750:                                VecMTDot_Seq,
751:                                VecScale_Seq,
752:                                VecCopy_Seq, /* 10 */
753:                                VecSet_Seq,
754:                                VecSwap_Seq,
755:                                VecAXPY_Seq,
756:                                VecAXPBY_Seq,
757:                                VecMAXPY_Seq,
758:                                VecAYPX_Seq,
759:                                VecWAXPY_Seq,
760:                                VecAXPBYPCZ_Seq,
761:                                VecPointwiseMult_Seq,
762:                                VecPointwiseDivide_Seq,
763:                                VecSetValues_Seq, /* 20 */
764:                                NULL,
765:                                NULL,
766:                                NULL,
767:                                VecGetSize_Seq,
768:                                VecGetSize_Seq,
769:                                NULL,
770:                                VecMax_Seq,
771:                                VecMin_Seq,
772:                                VecSetRandom_Seq,
773:                                VecSetOption_Seq, /* 30 */
774:                                VecSetValuesBlocked_Seq,
775:                                VecDestroy_Seq,
776:                                VecView_Seq,
777:                                VecPlaceArray_Seq,
778:                                VecReplaceArray_Seq,
779:                                VecDot_Seq,
780:                                VecTDot_Seq,
781:                                VecNorm_Seq,
782:                                VecMDot_Seq,
783:                                VecMTDot_Seq, /* 40 */
784:                                VecLoad_Default,
785:                                VecReciprocal_Default,
786:                                VecConjugate_Seq,
787:                                NULL,
788:                                NULL,
789:                                VecResetArray_Seq,
790:                                NULL,
791:                                VecMaxPointwiseDivide_Seq,
792:                                VecPointwiseMax_Seq,
793:                                VecPointwiseMaxAbs_Seq,
794:                                VecPointwiseMin_Seq,
795:                                VecGetValues_Seq,
796:                                NULL,
797:                                NULL,
798:                                NULL,
799:                                NULL,
800:                                NULL,
801:                                NULL,
802:                                VecStrideGather_Default,
803:                                VecStrideScatter_Default,
804:                                NULL,
805:                                NULL,
806:                                NULL,
807:                                NULL,
808:                                NULL,
809:                                VecStrideSubSetGather_Default,
810:                                VecStrideSubSetScatter_Default,
811:                                NULL,
812:                                NULL,
813:                                NULL
814: };


817: /*
818:       This is called by VecCreate_Seq() (i.e. VecCreateSeq()) and VecCreateSeqWithArray()
819: */
820: PetscErrorCode VecCreate_Seq_Private(Vec v,const PetscScalar array[])
821: {
822:   Vec_Seq        *s;

826:   PetscNewLog(v,&s);
827:   PetscMemcpy(v->ops,&DvOps,sizeof(DvOps));

829:   v->data            = (void*)s;
830:   v->petscnative     = PETSC_TRUE;
831:   s->array           = (PetscScalar*)array;
832:   s->array_allocated = NULL;
833:   if (array) v->offloadmask = PETSC_OFFLOAD_CPU;

835:   PetscLayoutSetUp(v->map);
836:   PetscObjectChangeTypeName((PetscObject)v,VECSEQ);
837: #if defined(PETSC_HAVE_MATLAB_ENGINE)
838:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEnginePut_C",VecMatlabEnginePut_Default);
839:   PetscObjectComposeFunction((PetscObject)v,"PetscMatlabEngineGet_C",VecMatlabEngineGet_Default);
840: #endif
841:   return(0);
842: }

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

848:    Collective

850:    Input Parameter:
851: +  comm - the communicator, should be PETSC_COMM_SELF
852: .  bs - the block size
853: .  n - the vector length
854: -  array - memory where the vector elements are to be stored.

856:    Output Parameter:
857: .  V - the vector

859:    Notes:
860:    Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
861:    same type as an existing vector.

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

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

869:    Level: intermediate

871: .seealso: VecCreateMPIWithArray(), VecCreate(), VecDuplicate(), VecDuplicateVecs(),
872:           VecCreateGhost(), VecCreateSeq(), VecPlaceArray()
873: @*/
874: PetscErrorCode  VecCreateSeqWithArray(MPI_Comm comm,PetscInt bs,PetscInt n,const PetscScalar array[],Vec *V)
875: {
877:   PetscMPIInt    size;

880:   VecCreate(comm,V);
881:   VecSetSizes(*V,n,n);
882:   VecSetBlockSize(*V,bs);
883:   MPI_Comm_size(comm,&size);
884:   if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than one process");
885:   VecCreate_Seq_Private(*V,array);
886:   return(0);
887: }