← Back to branch summary

~chaffra/+junk/petsc

« back to all changes in this revision

Viewing changes to .pc/upstream-p8.patch/src/vec/vec/utils/vinv.c

  • Committer: Package Import Robot
  • Author(s): "Adam C. Powell, IV", Adam C. Powell, IV, Matthias Klose
  • Date: 2011-04-07 21:20:50 UTC
  • mfrom: (6.1.17 sid)
  • Revision ID: package-import@ubuntu.com-20110407212050-gu6b1jbtvs8992nb
Tags: 3.1.dfsg-11
http://bugs.debian.org/602660
http://bugs.debian.org/605294
http://bugs.debian.org/621609
http://bugs.debian.org/608902
[ Adam C. Powell, IV ]
* Use HDF5 only on openmpi architectures (closes: #602660, #605294).
* Bump Standards-Version, no changes needed.
* Update to upstream patch level 8 (closes: #621609).
* Clean target update.
* Fix stamps to prevent double-building.

[ Matthias Klose ]
* Patch PETSc build system to link properly to hypre libs, and change BLACS
  lib order in rules, to make binutils-gold work (closes: #608902).
* Explicitly set PETSC_BUILD_GNU_SYSTEM in rules.

Show diffs side-by-side

added added

removed removed

Lines of Context:
.pc/upstream-p8.patch/src/vec/vec/utils/vinv.c
 
1
#define PETSCVEC_DLL
 
2
 
 
3
/*
 
4
     Some useful vector utility functions.
 
5
*/
 
6
#include "private/vecimpl.h"          /*I "petscvec.h" I*/
 
7
 
 
8
extern MPI_Op VecMax_Local_Op;
 
9
extern MPI_Op VecMin_Local_Op;
 
10
 
 
11
#undef __FUNCT__  
 
12
#define __FUNCT__ "VecStrideScale"
 
13
/*@
 
14
   VecStrideScale - Scales a subvector of a vector defined 
 
15
   by a starting point and a stride.
 
16
 
 
17
   Collective on Vec
 
18
 
 
19
   Input Parameter:
 
20
+  v - the vector 
 
21
.  start - starting point of the subvector (defined by a stride)
 
22
-  scale - value to multiply each subvector entry by
 
23
 
 
24
   Notes:
 
25
   One must call VecSetBlockSize() before this routine to set the stride 
 
26
   information, or use a vector created from a multicomponent DA.
 
27
 
 
28
   This will only work if the desire subvector is a stride subvector
 
29
 
 
30
   Level: advanced
 
31
 
 
32
   Concepts: scale^on stride of vector
 
33
   Concepts: stride^scale
 
34
 
 
35
.seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
 
36
@*/
 
37
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
 
38
{
 
39
  PetscErrorCode ierr;
 
40
  PetscInt       i,n,bs;
 
41
  PetscScalar    *x;
 
42
 
 
43
  PetscFunctionBegin;
 
44
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
45
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
46
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
47
 
 
48
  bs   = v->map->bs;
 
49
  if (start < 0) {
 
50
    SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 
51
  } else if (start >= bs) {
 
52
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n\
 
53
            Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
 
54
  }
 
55
  x += start;
 
56
 
 
57
  for (i=0; i<n; i+=bs) {
 
58
    x[i] *= scale;
 
59
  }
 
60
  x -= start;
 
61
 
 
62
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
63
  PetscFunctionReturn(0);
 
64
}
 
65
 
 
66
#undef __FUNCT__  
 
67
#define __FUNCT__ "VecStrideNorm"
 
68
/*@
 
69
   VecStrideNorm - Computes the norm of subvector of a vector defined 
 
70
   by a starting point and a stride.
 
71
 
 
72
   Collective on Vec
 
73
 
 
74
   Input Parameter:
 
75
+  v - the vector 
 
76
.  start - starting point of the subvector (defined by a stride)
 
77
-  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
 
78
 
 
79
   Output Parameter:
 
80
.  norm - the norm
 
81
 
 
82
   Notes:
 
83
   One must call VecSetBlockSize() before this routine to set the stride 
 
84
   information, or use a vector created from a multicomponent DA.
 
85
 
 
86
   If x is the array representing the vector x then this computes the norm 
 
87
   of the array (x[start],x[start+stride],x[start+2*stride], ....)
 
88
 
 
89
   This is useful for computing, say the norm of the pressure variable when
 
90
   the pressure is stored (interlaced) with other variables, say density etc.
 
91
 
 
92
   This will only work if the desire subvector is a stride subvector
 
93
 
 
94
   Level: advanced
 
95
 
 
96
   Concepts: norm^on stride of vector
 
97
   Concepts: stride^norm
 
98
 
 
99
.seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
 
100
@*/
 
101
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
 
102
{
 
103
  PetscErrorCode ierr;
 
104
  PetscInt       i,n,bs;
 
105
  PetscScalar    *x;
 
106
  PetscReal      tnorm;
 
107
  MPI_Comm       comm;
 
108
 
 
109
  PetscFunctionBegin;
 
110
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
111
  PetscValidDoublePointer(nrm,3);
 
112
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
113
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
114
  ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
 
115
 
 
116
  bs   = v->map->bs;
 
117
  if (start < 0) {
 
118
    SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 
119
  } else if (start >= bs) {
 
120
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n\
 
121
            Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
 
122
  }
 
123
  x += start;
 
124
 
 
125
  if (ntype == NORM_2) {
 
126
    PetscScalar sum = 0.0;
 
127
    for (i=0; i<n; i+=bs) {
 
128
      sum += x[i]*(PetscConj(x[i]));
 
129
    }
 
130
    tnorm  = PetscRealPart(sum);
 
131
    ierr   = MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_SUM,comm);CHKERRQ(ierr);
 
132
    *nrm = sqrt(*nrm);
 
133
  } else if (ntype == NORM_1) {
 
134
    tnorm = 0.0;
 
135
    for (i=0; i<n; i+=bs) {
 
136
      tnorm += PetscAbsScalar(x[i]);
 
137
    }
 
138
    ierr   = MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_SUM,comm);CHKERRQ(ierr);
 
139
  } else if (ntype == NORM_INFINITY) {
 
140
    PetscReal tmp;
 
141
    tnorm = 0.0;
 
142
 
 
143
    for (i=0; i<n; i+=bs) {
 
144
      if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
 
145
      /* check special case of tmp == NaN */
 
146
      if (tmp != tmp) {tnorm = tmp; break;}
 
147
    } 
 
148
    ierr   = MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_MAX,comm);CHKERRQ(ierr);
 
149
  } else {
 
150
    SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
 
151
  }
 
152
 
 
153
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
154
  PetscFunctionReturn(0);
 
155
}
 
156
 
 
157
#undef __FUNCT__  
 
158
#define __FUNCT__ "VecStrideMax"
 
159
/*@
 
160
   VecStrideMax - Computes the maximum of subvector of a vector defined 
 
161
   by a starting point and a stride and optionally its location.
 
162
 
 
163
   Collective on Vec
 
164
 
 
165
   Input Parameter:
 
166
+  v - the vector 
 
167
-  start - starting point of the subvector (defined by a stride)
 
168
 
 
169
   Output Parameter:
 
170
+  index - the location where the maximum occurred  (pass PETSC_NULL if not required)
 
171
-  nrm - the max
 
172
 
 
173
   Notes:
 
174
   One must call VecSetBlockSize() before this routine to set the stride 
 
175
   information, or use a vector created from a multicomponent DA.
 
176
 
 
177
   If xa is the array representing the vector x, then this computes the max
 
178
   of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
 
179
 
 
180
   This is useful for computing, say the maximum of the pressure variable when
 
181
   the pressure is stored (interlaced) with other variables, e.g., density, etc.
 
182
   This will only work if the desire subvector is a stride subvector.
 
183
 
 
184
   Level: advanced
 
185
 
 
186
   Concepts: maximum^on stride of vector
 
187
   Concepts: stride^maximum
 
188
 
 
189
.seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
 
190
@*/
 
191
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
 
192
{
 
193
  PetscErrorCode ierr;
 
194
  PetscInt       i,n,bs,id;
 
195
  PetscScalar    *x;
 
196
  PetscReal      max,tmp;
 
197
  MPI_Comm       comm;
 
198
 
 
199
  PetscFunctionBegin;
 
200
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
201
  PetscValidDoublePointer(nrm,3);
 
202
 
 
203
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
204
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
205
  ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
 
206
 
 
207
  bs   = v->map->bs;
 
208
  if (start < 0) {
 
209
    SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 
210
  } else if (start >= bs) {
 
211
    SETERRQ2(PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n\
 
212
            Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
 
213
  }
 
214
  x += start;
 
215
 
 
216
  id = -1;
 
217
  if (!n) {
 
218
    max = PETSC_MIN;
 
219
  } else {
 
220
    id  = 0;
 
221
#if defined(PETSC_USE_COMPLEX)
 
222
    max = PetscRealPart(x[0]);
 
223
#else
 
224
    max = x[0];
 
225
#endif
 
226
    for (i=bs; i<n; i+=bs) {
 
227
#if defined(PETSC_USE_COMPLEX)
 
228
      if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
 
229
#else
 
230
      if ((tmp = x[i]) > max) { max = tmp; id = i;} 
 
231
#endif
 
232
    }
 
233
  }
 
234
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
235
 
 
236
  if (!idex) {
 
237
    ierr   = MPI_Allreduce(&max,nrm,1,MPIU_REAL,MPI_MAX,comm);CHKERRQ(ierr);
 
238
  } else {
 
239
    PetscReal in[2],out[2];
 
240
    PetscInt  rstart;
 
241
 
 
242
    ierr  = VecGetOwnershipRange(v,&rstart,PETSC_NULL);CHKERRQ(ierr);
 
243
    in[0] = max;
 
244
    in[1] = rstart+id+start;
 
245
    ierr  = MPI_Allreduce(in,out,2,MPIU_REAL,VecMax_Local_Op,((PetscObject)v)->comm);CHKERRQ(ierr);
 
246
    *nrm  = out[0];
 
247
    *idex = (PetscInt)out[1];
 
248
  }
 
249
 
 
250
  PetscFunctionReturn(0);
 
251
}
 
252
 
 
253
#undef __FUNCT__  
 
254
#define __FUNCT__ "VecStrideMin"
 
255
/*@
 
256
   VecStrideMin - Computes the minimum of subvector of a vector defined 
 
257
   by a starting point and a stride and optionally its location.
 
258
 
 
259
   Collective on Vec
 
260
 
 
261
   Input Parameter:
 
262
+  v - the vector 
 
263
-  start - starting point of the subvector (defined by a stride)
 
264
 
 
265
   Output Parameter:
 
266
+  idex - the location where the minimum occurred. (pass PETSC_NULL if not required)
 
267
-  nrm - the min
 
268
 
 
269
   Level: advanced
 
270
 
 
271
   Notes:
 
272
   One must call VecSetBlockSize() before this routine to set the stride 
 
273
   information, or use a vector created from a multicomponent DA.
 
274
 
 
275
   If xa is the array representing the vector x, then this computes the min
 
276
   of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)
 
277
 
 
278
   This is useful for computing, say the minimum of the pressure variable when
 
279
   the pressure is stored (interlaced) with other variables, e.g., density, etc.
 
280
   This will only work if the desire subvector is a stride subvector.
 
281
 
 
282
   Concepts: minimum^on stride of vector
 
283
   Concepts: stride^minimum
 
284
 
 
285
.seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
 
286
@*/
 
287
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
 
288
{
 
289
  PetscErrorCode ierr;
 
290
  PetscInt       i,n,bs,id;
 
291
  PetscScalar    *x;
 
292
  PetscReal      min,tmp;
 
293
  MPI_Comm       comm;
 
294
 
 
295
  PetscFunctionBegin;
 
296
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
297
  PetscValidDoublePointer(nrm,4);
 
298
 
 
299
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
300
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
301
  ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
 
302
 
 
303
  bs   = v->map->bs;
 
304
  if (start < 0) {
 
305
    SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 
306
  } else if (start >= bs) {
 
307
    SETERRQ2(PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n\
 
308
            Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
 
309
  }
 
310
  x += start;
 
311
 
 
312
  id = -1;
 
313
  if (!n) {
 
314
    min = PETSC_MAX;
 
315
  } else {
 
316
    id = 0;
 
317
#if defined(PETSC_USE_COMPLEX)
 
318
    min = PetscRealPart(x[0]);
 
319
#else
 
320
    min = x[0];
 
321
#endif
 
322
    for (i=bs; i<n; i+=bs) {
 
323
#if defined(PETSC_USE_COMPLEX)
 
324
      if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
 
325
#else
 
326
      if ((tmp = x[i]) < min) { min = tmp; id = i;} 
 
327
#endif
 
328
    }
 
329
  }
 
330
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
331
 
 
332
  if (!idex) {
 
333
    ierr   = MPI_Allreduce(&min,nrm,1,MPIU_REAL,MPI_MIN,comm);CHKERRQ(ierr);
 
334
  } else {
 
335
    PetscReal in[2],out[2];
 
336
    PetscInt  rstart;
 
337
 
 
338
    ierr  = VecGetOwnershipRange(v,&rstart,PETSC_NULL);CHKERRQ(ierr);
 
339
    in[0] = min;
 
340
    in[1] = rstart+id;
 
341
    ierr  = MPI_Allreduce(in,out,2,MPIU_REAL,VecMin_Local_Op,((PetscObject)v)->comm);CHKERRQ(ierr);
 
342
    *nrm  = out[0];
 
343
    *idex = (PetscInt)out[1];
 
344
  }
 
345
 
 
346
  PetscFunctionReturn(0);
 
347
}
 
348
 
 
349
#undef __FUNCT__  
 
350
#define __FUNCT__ "VecStrideScaleAll"
 
351
/*@
 
352
   VecStrideScaleAll - Scales the subvectors of a vector defined 
 
353
   by a starting point and a stride.
 
354
 
 
355
   Collective on Vec
 
356
 
 
357
   Input Parameter:
 
358
+  v - the vector 
 
359
-  scales - values to multiply each subvector entry by
 
360
 
 
361
   Notes:
 
362
   One must call VecSetBlockSize() before this routine to set the stride 
 
363
   information, or use a vector created from a multicomponent DA.
 
364
 
 
365
 
 
366
   Level: advanced
 
367
 
 
368
   Concepts: scale^on stride of vector
 
369
   Concepts: stride^scale
 
370
 
 
371
.seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
 
372
@*/
 
373
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideScaleAll(Vec v,PetscScalar *scales)
 
374
{
 
375
  PetscErrorCode ierr;
 
376
  PetscInt       i,j,n,bs;
 
377
  PetscScalar    *x;
 
378
 
 
379
  PetscFunctionBegin;
 
380
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
381
  PetscValidScalarPointer(scales,2);
 
382
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
383
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
384
 
 
385
  bs   = v->map->bs;
 
386
 
 
387
  /* need to provide optimized code for each bs */
 
388
  for (i=0; i<n; i+=bs) {
 
389
    for (j=0; j<bs; j++) {
 
390
      x[i+j] *= scales[j];
 
391
    }
 
392
  }
 
393
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
394
  PetscFunctionReturn(0);
 
395
}
 
396
 
 
397
#undef __FUNCT__  
 
398
#define __FUNCT__ "VecStrideNormAll"
 
399
/*@
 
400
   VecStrideNormAll - Computes the norms  subvectors of a vector defined 
 
401
   by a starting point and a stride.
 
402
 
 
403
   Collective on Vec
 
404
 
 
405
   Input Parameter:
 
406
+  v - the vector 
 
407
-  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY
 
408
 
 
409
   Output Parameter:
 
410
.  nrm - the norms
 
411
 
 
412
   Notes:
 
413
   One must call VecSetBlockSize() before this routine to set the stride 
 
414
   information, or use a vector created from a multicomponent DA.
 
415
 
 
416
   If x is the array representing the vector x then this computes the norm 
 
417
   of the array (x[start],x[start+stride],x[start+2*stride], ....)
 
418
 
 
419
   This is useful for computing, say the norm of the pressure variable when
 
420
   the pressure is stored (interlaced) with other variables, say density etc.
 
421
 
 
422
   This will only work if the desire subvector is a stride subvector
 
423
 
 
424
   Level: advanced
 
425
 
 
426
   Concepts: norm^on stride of vector
 
427
   Concepts: stride^norm
 
428
 
 
429
.seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
 
430
@*/
 
431
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
 
432
{
 
433
  PetscErrorCode ierr;
 
434
  PetscInt       i,j,n,bs;
 
435
  PetscScalar    *x;
 
436
  PetscReal      tnorm[128];
 
437
  MPI_Comm       comm;
 
438
 
 
439
  PetscFunctionBegin;
 
440
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
441
  PetscValidDoublePointer(nrm,3);
 
442
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
443
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
444
  ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
 
445
 
 
446
  bs   = v->map->bs;
 
447
  if (bs > 128) SETERRQ(PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
 
448
 
 
449
  if (ntype == NORM_2) {
 
450
    PetscScalar sum[128];
 
451
    for (j=0; j<bs; j++) sum[j] = 0.0;
 
452
    for (i=0; i<n; i+=bs) {
 
453
      for (j=0; j<bs; j++) {
 
454
        sum[j] += x[i+j]*(PetscConj(x[i+j]));
 
455
      }
 
456
    }
 
457
    for (j=0; j<bs; j++) {
 
458
      tnorm[j]  = PetscRealPart(sum[j]);
 
459
    }
 
460
    ierr   = MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPI_SUM,comm);CHKERRQ(ierr);
 
461
    for (j=0; j<bs; j++) {
 
462
      nrm[j] = sqrt(nrm[j]);
 
463
    }
 
464
  } else if (ntype == NORM_1) {
 
465
    for (j=0; j<bs; j++) {
 
466
      tnorm[j] = 0.0;
 
467
    }
 
468
    for (i=0; i<n; i+=bs) {
 
469
      for (j=0; j<bs; j++) {
 
470
        tnorm[j] += PetscAbsScalar(x[i+j]);
 
471
      }
 
472
    }
 
473
    ierr   = MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPI_SUM,comm);CHKERRQ(ierr);
 
474
  } else if (ntype == NORM_INFINITY) {
 
475
    PetscReal tmp;
 
476
    for (j=0; j<bs; j++) {
 
477
      tnorm[j] = 0.0;
 
478
    }
 
479
 
 
480
    for (i=0; i<n; i+=bs) {
 
481
      for (j=0; j<bs; j++) {
 
482
        if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
 
483
        /* check special case of tmp == NaN */
 
484
        if (tmp != tmp) {tnorm[j] = tmp; break;}
 
485
      }
 
486
    } 
 
487
    ierr   = MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPI_MAX,comm);CHKERRQ(ierr);
 
488
  } else {
 
489
    SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
 
490
  }
 
491
 
 
492
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
493
  PetscFunctionReturn(0);
 
494
}
 
495
 
 
496
#undef __FUNCT__  
 
497
#define __FUNCT__ "VecStrideMaxAll"
 
498
/*@
 
499
   VecStrideMaxAll - Computes the maximums of subvectors of a vector defined 
 
500
   by a starting point and a stride and optionally its location.
 
501
 
 
502
   Collective on Vec
 
503
 
 
504
   Input Parameter:
 
505
.  v - the vector 
 
506
 
 
507
   Output Parameter:
 
508
+  index - the location where the maximum occurred (not supported, pass PETSC_NULL,
 
509
           if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
 
510
-  nrm - the maximums
 
511
 
 
512
   Notes:
 
513
   One must call VecSetBlockSize() before this routine to set the stride 
 
514
   information, or use a vector created from a multicomponent DA.
 
515
 
 
516
   This is useful for computing, say the maximum of the pressure variable when
 
517
   the pressure is stored (interlaced) with other variables, e.g., density, etc.
 
518
   This will only work if the desire subvector is a stride subvector.
 
519
 
 
520
   Level: advanced
 
521
 
 
522
   Concepts: maximum^on stride of vector
 
523
   Concepts: stride^maximum
 
524
 
 
525
.seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
 
526
@*/
 
527
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
 
528
{
 
529
  PetscErrorCode ierr;
 
530
  PetscInt       i,j,n,bs;
 
531
  PetscScalar    *x;
 
532
  PetscReal      max[128],tmp;
 
533
  MPI_Comm       comm;
 
534
 
 
535
  PetscFunctionBegin;
 
536
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
537
  PetscValidDoublePointer(nrm,3);
 
538
  if (idex) {
 
539
    SETERRQ(PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
 
540
  }
 
541
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
542
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
543
  ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
 
544
 
 
545
  bs   = v->map->bs;
 
546
  if (bs > 128) SETERRQ(PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
 
547
 
 
548
  if (!n) {
 
549
    for (j=0; j<bs; j++) {
 
550
      max[j] = PETSC_MIN;
 
551
    }
 
552
  } else {
 
553
    for (j=0; j<bs; j++) {
 
554
#if defined(PETSC_USE_COMPLEX)
 
555
      max[j] = PetscRealPart(x[j]);
 
556
#else
 
557
      max[j] = x[j];
 
558
#endif
 
559
    }
 
560
    for (i=bs; i<n; i+=bs) {
 
561
      for (j=0; j<bs; j++) {
 
562
#if defined(PETSC_USE_COMPLEX)
 
563
        if ((tmp = PetscRealPart(x[i+j])) > max[j]) { max[j] = tmp;}
 
564
#else
 
565
        if ((tmp = x[i+j]) > max[j]) { max[j] = tmp; } 
 
566
#endif
 
567
      }
 
568
    }
 
569
  }
 
570
  ierr   = MPI_Allreduce(max,nrm,bs,MPIU_REAL,MPI_MAX,comm);CHKERRQ(ierr);
 
571
 
 
572
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
573
  PetscFunctionReturn(0);
 
574
}
 
575
 
 
576
#undef __FUNCT__  
 
577
#define __FUNCT__ "VecStrideMinAll"
 
578
/*@
 
579
   VecStrideMinAll - Computes the minimum of subvector of a vector defined 
 
580
   by a starting point and a stride and optionally its location.
 
581
 
 
582
   Collective on Vec
 
583
 
 
584
   Input Parameter:
 
585
.  v - the vector 
 
586
 
 
587
   Output Parameter:
 
588
+  idex - the location where the minimum occurred (not supported, pass PETSC_NULL,
 
589
           if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
 
590
-  nrm - the minimums
 
591
 
 
592
   Level: advanced
 
593
 
 
594
   Notes:
 
595
   One must call VecSetBlockSize() before this routine to set the stride 
 
596
   information, or use a vector created from a multicomponent DA.
 
597
 
 
598
   This is useful for computing, say the minimum of the pressure variable when
 
599
   the pressure is stored (interlaced) with other variables, e.g., density, etc.
 
600
   This will only work if the desire subvector is a stride subvector.
 
601
 
 
602
   Concepts: minimum^on stride of vector
 
603
   Concepts: stride^minimum
 
604
 
 
605
.seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
 
606
@*/
 
607
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
 
608
{
 
609
  PetscErrorCode ierr;
 
610
  PetscInt       i,n,bs,j;
 
611
  PetscScalar    *x;
 
612
  PetscReal      min[128],tmp;
 
613
  MPI_Comm       comm;
 
614
 
 
615
  PetscFunctionBegin;
 
616
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
617
  PetscValidDoublePointer(nrm,3);
 
618
  if (idex) {
 
619
    SETERRQ(PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
 
620
  }
 
621
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
622
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
623
  ierr = PetscObjectGetComm((PetscObject)v,&comm);CHKERRQ(ierr);
 
624
 
 
625
  bs   = v->map->bs;
 
626
  if (bs > 128) SETERRQ(PETSC_ERR_SUP,"Currently supports only blocksize up to 128");
 
627
 
 
628
  if (!n) {
 
629
    for (j=0; j<bs; j++) {
 
630
      min[j] = PETSC_MAX;
 
631
    }
 
632
  } else {
 
633
    for (j=0; j<bs; j++) {
 
634
#if defined(PETSC_USE_COMPLEX)
 
635
      min[j] = PetscRealPart(x[j]);
 
636
#else
 
637
      min[j] = x[j];
 
638
#endif
 
639
    }
 
640
    for (i=bs; i<n; i+=bs) {
 
641
      for (j=0; j<bs; j++) {
 
642
#if defined(PETSC_USE_COMPLEX)
 
643
        if ((tmp = PetscRealPart(x[i+j])) < min[j]) { min[j] = tmp;}
 
644
#else
 
645
        if ((tmp = x[i+j]) < min[j]) { min[j] = tmp; } 
 
646
#endif
 
647
      }
 
648
    }
 
649
  }
 
650
  ierr   = MPI_Allreduce(min,nrm,bs,MPIU_REAL,MPI_MIN,comm);CHKERRQ(ierr);
 
651
 
 
652
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
653
  PetscFunctionReturn(0);
 
654
}
 
655
 
 
656
/*----------------------------------------------------------------------------------------------*/
 
657
#undef __FUNCT__  
 
658
#define __FUNCT__ "VecStrideGatherAll"
 
659
/*@
 
660
   VecStrideGatherAll - Gathers all the single components from a multi-component vector into
 
661
   separate vectors.
 
662
 
 
663
   Collective on Vec
 
664
 
 
665
   Input Parameter:
 
666
+  v - the vector 
 
667
-  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
 
668
 
 
669
   Output Parameter:
 
670
.  s - the location where the subvectors are stored
 
671
 
 
672
   Notes:
 
673
   One must call VecSetBlockSize() before this routine to set the stride 
 
674
   information, or use a vector created from a multicomponent DA.
 
675
 
 
676
   If x is the array representing the vector x then this gathers
 
677
   the arrays (x[start],x[start+stride],x[start+2*stride], ....)
 
678
   for start=0,1,2,...bs-1
 
679
 
 
680
   The parallel layout of the vector and the subvector must be the same;
 
681
   i.e., nlocal of v = stride*(nlocal of s) 
 
682
 
 
683
   Not optimized; could be easily
 
684
 
 
685
   Level: advanced
 
686
 
 
687
   Concepts: gather^into strided vector
 
688
 
 
689
.seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
 
690
          VecStrideScatterAll()
 
691
@*/
 
692
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
 
693
{
 
694
  PetscErrorCode ierr;
 
695
  PetscInt       i,n,n2,bs,j,k,*bss = PETSC_NULL,nv,jj,nvc;
 
696
  PetscScalar    *x,**y;
 
697
 
 
698
  PetscFunctionBegin;
 
699
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
700
  PetscValidPointer(s,2);
 
701
  PetscValidHeaderSpecific(*s,VEC_COOKIE,2);
 
702
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
703
  ierr = VecGetLocalSize(s[0],&n2);CHKERRQ(ierr);
 
704
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
705
  bs   = v->map->bs;
 
706
  if (bs < 0) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
 
707
 
 
708
  ierr = PetscMalloc2(bs,PetscReal*,&y,bs,PetscInt,&bss);CHKERRQ(ierr);
 
709
  nv   = 0;
 
710
  nvc  = 0;
 
711
  for (i=0; i<bs; i++) {
 
712
    ierr = VecGetBlockSize(s[i],&bss[i]);CHKERRQ(ierr);
 
713
    if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
 
714
    ierr = VecGetArray(s[i],&y[i]);CHKERRQ(ierr);
 
715
    nvc  += bss[i];
 
716
    nv++;
 
717
    if (nvc > bs)  SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
 
718
    if (nvc == bs) break;
 
719
  }
 
720
 
 
721
  n =  n/bs;
 
722
 
 
723
  jj = 0;
 
724
  if (addv == INSERT_VALUES) {
 
725
    for (j=0; j<nv; j++) {
 
726
      for (k=0; k<bss[j]; k++) {
 
727
        for (i=0; i<n; i++) {
 
728
          y[j][i*bss[j] + k] = x[bs*i+jj+k];
 
729
        }
 
730
      }
 
731
      jj += bss[j];
 
732
    }
 
733
  } else if (addv == ADD_VALUES) {
 
734
    for (j=0; j<nv; j++) {
 
735
      for (k=0; k<bss[j]; k++) {
 
736
        for (i=0; i<n; i++) {
 
737
          y[j][i*bss[j] + k] += x[bs*i+jj+k];
 
738
        }
 
739
      }
 
740
      jj += bss[j];
 
741
    }
 
742
#if !defined(PETSC_USE_COMPLEX)
 
743
  } else if (addv == MAX_VALUES) {
 
744
    for (j=0; j<nv; j++) {
 
745
      for (k=0; k<bss[j]; k++) {
 
746
        for (i=0; i<n; i++) {
 
747
          y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
 
748
        }
 
749
      }
 
750
      jj += bss[j];
 
751
    }
 
752
#endif
 
753
  } else {
 
754
    SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
 
755
  }
 
756
 
 
757
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
758
  for (i=0; i<nv; i++) {
 
759
    ierr = VecRestoreArray(s[i],&y[i]);CHKERRQ(ierr);
 
760
  }
 
761
 
 
762
  ierr = PetscFree2(y,bss);CHKERRQ(ierr);
 
763
  PetscFunctionReturn(0);
 
764
}
 
765
 
 
766
#undef __FUNCT__  
 
767
#define __FUNCT__ "VecStrideScatterAll"
 
768
/*@
 
769
   VecStrideScatterAll - Scatters all the single components from separate vectors into 
 
770
     a multi-component vector.
 
771
 
 
772
   Collective on Vec
 
773
 
 
774
   Input Parameter:
 
775
+  s - the location where the subvectors are stored
 
776
-  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
 
777
 
 
778
   Output Parameter:
 
779
.  v - the multicomponent vector 
 
780
 
 
781
   Notes:
 
782
   One must call VecSetBlockSize() before this routine to set the stride 
 
783
   information, or use a vector created from a multicomponent DA.
 
784
 
 
785
   The parallel layout of the vector and the subvector must be the same;
 
786
   i.e., nlocal of v = stride*(nlocal of s) 
 
787
 
 
788
   Not optimized; could be easily
 
789
 
 
790
   Level: advanced
 
791
 
 
792
   Concepts:  scatter^into strided vector
 
793
 
 
794
.seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
 
795
          VecStrideScatterAll()
 
796
@*/
 
797
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
 
798
{
 
799
  PetscErrorCode ierr;
 
800
  PetscInt        i,n,n2,bs,j,jj,k,*bss = PETSC_NULL,nv,nvc;
 
801
  PetscScalar     *x,**y;
 
802
 
 
803
  PetscFunctionBegin;
 
804
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
805
  PetscValidPointer(s,2);
 
806
  PetscValidHeaderSpecific(*s,VEC_COOKIE,2);
 
807
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
808
  ierr = VecGetLocalSize(s[0],&n2);CHKERRQ(ierr);
 
809
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
810
  bs   = v->map->bs;
 
811
  if (bs < 0) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");
 
812
 
 
813
  ierr = PetscMalloc2(bs,PetscScalar**,&y,bs,PetscInt,&bss);CHKERRQ(ierr);
 
814
  nv  = 0;
 
815
  nvc = 0;
 
816
  for (i=0; i<bs; i++) {
 
817
    ierr = VecGetBlockSize(s[i],&bss[i]);CHKERRQ(ierr);
 
818
    if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
 
819
    ierr = VecGetArray(s[i],&y[i]);CHKERRQ(ierr);
 
820
    nvc  += bss[i];
 
821
    nv++;
 
822
    if (nvc > bs)  SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
 
823
    if (nvc == bs) break;
 
824
  }
 
825
 
 
826
  n =  n/bs;
 
827
 
 
828
  jj = 0;
 
829
  if (addv == INSERT_VALUES) {
 
830
    for (j=0; j<nv; j++) {
 
831
      for (k=0; k<bss[j]; k++) {
 
832
        for (i=0; i<n; i++) {
 
833
          x[bs*i+jj+k] = y[j][i*bss[j] + k];
 
834
        }
 
835
      }
 
836
      jj += bss[j];
 
837
    }
 
838
  } else if (addv == ADD_VALUES) {
 
839
    for (j=0; j<nv; j++) {
 
840
      for (k=0; k<bss[j]; k++) {
 
841
        for (i=0; i<n; i++) {
 
842
          x[bs*i+jj+k] += y[j][i*bss[j] + k];
 
843
        }
 
844
      }
 
845
      jj += bss[j];
 
846
    }
 
847
#if !defined(PETSC_USE_COMPLEX)
 
848
  } else if (addv == MAX_VALUES) {
 
849
    for (j=0; j<nv; j++) {
 
850
      for (k=0; k<bss[j]; k++) {
 
851
        for (i=0; i<n; i++) {
 
852
          x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
 
853
        }
 
854
      }
 
855
      jj += bss[j];
 
856
    }
 
857
#endif
 
858
  } else {
 
859
    SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
 
860
  }
 
861
 
 
862
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
863
  for (i=0; i<nv; i++) {
 
864
    ierr = VecRestoreArray(s[i],&y[i]);CHKERRQ(ierr);
 
865
  }
 
866
  ierr = PetscFree2(y,bss);CHKERRQ(ierr);
 
867
  PetscFunctionReturn(0);
 
868
}
 
869
 
 
870
#undef __FUNCT__  
 
871
#define __FUNCT__ "VecStrideGather"
 
872
/*@
 
873
   VecStrideGather - Gathers a single component from a multi-component vector into
 
874
   another vector.
 
875
 
 
876
   Collective on Vec
 
877
 
 
878
   Input Parameter:
 
879
+  v - the vector 
 
880
.  start - starting point of the subvector (defined by a stride)
 
881
-  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
 
882
 
 
883
   Output Parameter:
 
884
.  s - the location where the subvector is stored
 
885
 
 
886
   Notes:
 
887
   One must call VecSetBlockSize() before this routine to set the stride 
 
888
   information, or use a vector created from a multicomponent DA.
 
889
 
 
890
   If x is the array representing the vector x then this gathers
 
891
   the array (x[start],x[start+stride],x[start+2*stride], ....)
 
892
 
 
893
   The parallel layout of the vector and the subvector must be the same;
 
894
   i.e., nlocal of v = stride*(nlocal of s) 
 
895
 
 
896
   Not optimized; could be easily
 
897
 
 
898
   Level: advanced
 
899
 
 
900
   Concepts: gather^into strided vector
 
901
 
 
902
.seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
 
903
          VecStrideScatterAll()
 
904
@*/
 
905
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
 
906
{
 
907
  PetscErrorCode ierr;
 
908
  PetscInt       i,n,bs,ns;
 
909
  PetscScalar    *x,*y;
 
910
 
 
911
  PetscFunctionBegin;
 
912
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
913
  PetscValidHeaderSpecific(s,VEC_COOKIE,3);
 
914
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
915
  ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
 
916
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
917
  ierr = VecGetArray(s,&y);CHKERRQ(ierr);
 
918
 
 
919
  bs   = v->map->bs;
 
920
  if (start < 0) {
 
921
    SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 
922
  } else if (start >= bs) {
 
923
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n\
 
924
            Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
 
925
  }
 
926
  if (n != ns*bs) {
 
927
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for gather from original vector %D",ns*bs,n);
 
928
  }
 
929
  x += start;
 
930
  n =  n/bs;
 
931
 
 
932
  if (addv == INSERT_VALUES) {
 
933
    for (i=0; i<n; i++) {
 
934
      y[i] = x[bs*i];
 
935
    }
 
936
  } else if (addv == ADD_VALUES) {
 
937
    for (i=0; i<n; i++) {
 
938
      y[i] += x[bs*i];
 
939
    }
 
940
#if !defined(PETSC_USE_COMPLEX)
 
941
  } else if (addv == MAX_VALUES) {
 
942
    for (i=0; i<n; i++) {
 
943
      y[i] = PetscMax(y[i],x[bs*i]);
 
944
    }
 
945
#endif
 
946
  } else {
 
947
    SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
 
948
  }
 
949
 
 
950
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
951
  ierr = VecRestoreArray(s,&y);CHKERRQ(ierr);
 
952
  PetscFunctionReturn(0);
 
953
}
 
954
 
 
955
#undef __FUNCT__  
 
956
#define __FUNCT__ "VecStrideScatter"
 
957
/*@
 
958
   VecStrideScatter - Scatters a single component from a vector into a multi-component vector.
 
959
 
 
960
   Collective on Vec
 
961
 
 
962
   Input Parameter:
 
963
+  s - the single-component vector 
 
964
.  start - starting point of the subvector (defined by a stride)
 
965
-  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES
 
966
 
 
967
   Output Parameter:
 
968
.  v - the location where the subvector is scattered (the multi-component vector)
 
969
 
 
970
   Notes:
 
971
   One must call VecSetBlockSize() on the multi-component vector before this
 
972
   routine to set the stride  information, or use a vector created from a multicomponent DA.
 
973
 
 
974
   The parallel layout of the vector and the subvector must be the same;
 
975
   i.e., nlocal of v = stride*(nlocal of s) 
 
976
 
 
977
   Not optimized; could be easily
 
978
 
 
979
   Level: advanced
 
980
 
 
981
   Concepts: scatter^into strided vector
 
982
 
 
983
.seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
 
984
          VecStrideScatterAll()
 
985
@*/
 
986
PetscErrorCode PETSCVEC_DLLEXPORT VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
 
987
{
 
988
  PetscErrorCode ierr;
 
989
  PetscInt       i,n,bs,ns;
 
990
  PetscScalar    *x,*y;
 
991
 
 
992
  PetscFunctionBegin;
 
993
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
994
  PetscValidHeaderSpecific(s,VEC_COOKIE,3);
 
995
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
996
  ierr = VecGetLocalSize(s,&ns);CHKERRQ(ierr);
 
997
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
998
  ierr = VecGetArray(s,&y);CHKERRQ(ierr);
 
999
 
 
1000
  bs   = v->map->bs;
 
1001
  if (start < 0) {
 
1002
    SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
 
1003
  } else if (start >= bs) {
 
1004
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n\
 
1005
            Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
 
1006
  }
 
1007
  if (n != ns*bs) {
 
1008
    SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for scatter to multicomponent vector %D",ns*bs,n);
 
1009
  }
 
1010
  x += start;
 
1011
  n =  n/bs;
 
1012
 
 
1013
 
 
1014
  if (addv == INSERT_VALUES) {
 
1015
    for (i=0; i<n; i++) {
 
1016
      x[bs*i] = y[i];
 
1017
    }
 
1018
  } else if (addv == ADD_VALUES) {
 
1019
    for (i=0; i<n; i++) {
 
1020
      x[bs*i] += y[i];
 
1021
    }
 
1022
#if !defined(PETSC_USE_COMPLEX)
 
1023
  } else if (addv == MAX_VALUES) {
 
1024
    for (i=0; i<n; i++) {
 
1025
      x[bs*i] = PetscMax(y[i],x[bs*i]);
 
1026
    }
 
1027
#endif
 
1028
  } else {
 
1029
    SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
 
1030
  }
 
1031
 
 
1032
 
 
1033
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
1034
  ierr = VecRestoreArray(s,&y);CHKERRQ(ierr);
 
1035
  PetscFunctionReturn(0);
 
1036
}
 
1037
 
 
1038
#undef __FUNCT__  
 
1039
#define __FUNCT__ "VecReciprocal_Default"
 
1040
PetscErrorCode VecReciprocal_Default(Vec v)
 
1041
{
 
1042
  PetscErrorCode ierr;
 
1043
  PetscInt       i,n;
 
1044
  PetscScalar    *x;
 
1045
 
 
1046
  PetscFunctionBegin;
 
1047
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
1048
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
1049
  for (i=0; i<n; i++) {
 
1050
    if (x[i] != 0.0) x[i] = 1.0/x[i];
 
1051
  }
 
1052
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
1053
  PetscFunctionReturn(0);
 
1054
}
 
1055
 
 
1056
#undef __FUNCT__
 
1057
#define __FUNCT__ "VecExp"
 
1058
/*@
 
1059
  VecExp - Replaces each component of a vector by e^x_i
 
1060
 
 
1061
  Not collective
 
1062
 
 
1063
  Input Parameter:
 
1064
. v - The vector
 
1065
 
 
1066
  Output Parameter:
 
1067
. v - The vector of exponents
 
1068
 
 
1069
  Level: beginner
 
1070
 
 
1071
.seealso:  VecLog(), VecAbs(), VecSqrt(), VecReciprocal()
 
1072
 
 
1073
.keywords: vector, sqrt, square root
 
1074
@*/
 
1075
PetscErrorCode PETSCVEC_DLLEXPORT VecExp(Vec v)
 
1076
{
 
1077
  PetscScalar    *x;
 
1078
  PetscInt       i, n;
 
1079
  PetscErrorCode ierr;
 
1080
 
 
1081
  PetscFunctionBegin;
 
1082
  PetscValidHeaderSpecific(v, VEC_COOKIE,1);
 
1083
  if (v->ops->exp) {
 
1084
    ierr = (*v->ops->exp)(v);CHKERRQ(ierr);
 
1085
  } else {
 
1086
    ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
 
1087
    ierr = VecGetArray(v, &x);CHKERRQ(ierr);
 
1088
    for(i = 0; i < n; i++) {
 
1089
      x[i] = PetscExpScalar(x[i]);
 
1090
    }
 
1091
    ierr = VecRestoreArray(v, &x);CHKERRQ(ierr);
 
1092
  }
 
1093
  PetscFunctionReturn(0);
 
1094
}
 
1095
 
 
1096
#undef __FUNCT__
 
1097
#define __FUNCT__ "VecLog"
 
1098
/*@
 
1099
  VecLog - Replaces each component of a vector by log(x_i), the natural logarithm
 
1100
 
 
1101
  Not collective
 
1102
 
 
1103
  Input Parameter:
 
1104
. v - The vector
 
1105
 
 
1106
  Output Parameter:
 
1107
. v - The vector of logs
 
1108
 
 
1109
  Level: beginner
 
1110
 
 
1111
.seealso:  VecExp(), VecAbs(), VecSqrt(), VecReciprocal()
 
1112
 
 
1113
.keywords: vector, sqrt, square root
 
1114
@*/
 
1115
PetscErrorCode PETSCVEC_DLLEXPORT VecLog(Vec v)
 
1116
{
 
1117
  PetscScalar    *x;
 
1118
  PetscInt       i, n;
 
1119
  PetscErrorCode ierr;
 
1120
 
 
1121
  PetscFunctionBegin;
 
1122
  PetscValidHeaderSpecific(v, VEC_COOKIE,1);
 
1123
  if (v->ops->log) {
 
1124
    ierr = (*v->ops->log)(v);CHKERRQ(ierr);
 
1125
  } else {
 
1126
    ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
 
1127
    ierr = VecGetArray(v, &x);CHKERRQ(ierr);
 
1128
    for(i = 0; i < n; i++) {
 
1129
      x[i] = PetscLogScalar(x[i]);
 
1130
    }
 
1131
    ierr = VecRestoreArray(v, &x);CHKERRQ(ierr);
 
1132
  }
 
1133
  PetscFunctionReturn(0);
 
1134
}
 
1135
 
 
1136
#undef __FUNCT__
 
1137
#define __FUNCT__ "VecSqrt"
 
1138
/*@
 
1139
  VecSqrt - Replaces each component of a vector by the square root of its magnitude.
 
1140
 
 
1141
  Not collective
 
1142
 
 
1143
  Input Parameter:
 
1144
. v - The vector
 
1145
 
 
1146
  Output Parameter:
 
1147
. v - The vector square root
 
1148
 
 
1149
  Level: beginner
 
1150
 
 
1151
  Note: The actual function is sqrt(|x_i|)
 
1152
 
 
1153
.seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()
 
1154
 
 
1155
.keywords: vector, sqrt, square root
 
1156
@*/
 
1157
PetscErrorCode PETSCVEC_DLLEXPORT VecSqrt(Vec v)
 
1158
{
 
1159
  PetscScalar    *x;
 
1160
  PetscInt       i, n;
 
1161
  PetscErrorCode ierr;
 
1162
 
 
1163
  PetscFunctionBegin;
 
1164
  PetscValidHeaderSpecific(v, VEC_COOKIE,1);
 
1165
  if (v->ops->sqrt) {
 
1166
    ierr = (*v->ops->sqrt)(v);CHKERRQ(ierr);
 
1167
  } else {
 
1168
    ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
 
1169
    ierr = VecGetArray(v, &x);CHKERRQ(ierr);
 
1170
    for(i = 0; i < n; i++) {
 
1171
      x[i] = sqrt(PetscAbsScalar(x[i]));
 
1172
    }
 
1173
    ierr = VecRestoreArray(v, &x);CHKERRQ(ierr);
 
1174
  }
 
1175
  PetscFunctionReturn(0);
 
1176
}
 
1177
 
 
1178
#undef __FUNCT__
 
1179
#define __FUNCT__ "VecDotNorm2"
 
1180
/*@
 
1181
  VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector
 
1182
 
 
1183
  Collective on Vec
 
1184
 
 
1185
  Input Parameter:
 
1186
+ s - first vector
 
1187
- t - second vector
 
1188
 
 
1189
  Output Parameter:
 
1190
+ dp - s't
 
1191
- nm - t't
 
1192
 
 
1193
  Level: advanced
 
1194
 
 
1195
.seealso:   VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()
 
1196
 
 
1197
.keywords: vector, sqrt, square root
 
1198
@*/
 
1199
PetscErrorCode PETSCVEC_DLLEXPORT VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscScalar *nm)
 
1200
{
 
1201
  PetscScalar    *sx, *tx, dpx = 0.0, nmx = 0.0,work[2],sum[2];
 
1202
  PetscInt       i, n;
 
1203
  PetscErrorCode ierr;
 
1204
 
 
1205
  PetscFunctionBegin;
 
1206
  PetscValidHeaderSpecific(s, VEC_COOKIE,1);
 
1207
  PetscValidHeaderSpecific(t, VEC_COOKIE,2);
 
1208
  PetscValidScalarPointer(dp,3);
 
1209
  PetscValidScalarPointer(nm,4);
 
1210
  PetscValidType(s,1);
 
1211
  PetscValidType(t,2);
 
1212
  PetscCheckSameTypeAndComm(s,1,t,2);
 
1213
  if (s->map->N != t->map->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
 
1214
  if (s->map->n != t->map->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
 
1215
 
 
1216
  ierr = PetscLogEventBarrierBegin(VEC_DotNormBarrier,s,t,0,0,((PetscObject)s)->comm);CHKERRQ(ierr);
 
1217
  if (s->ops->dotnorm2) {
 
1218
    ierr = (*s->ops->dotnorm2)(s,t,dp,nm);CHKERRQ(ierr);
 
1219
  } else {
 
1220
    ierr = VecGetLocalSize(s, &n);CHKERRQ(ierr);
 
1221
    ierr = VecGetArray(s, &sx);CHKERRQ(ierr);
 
1222
    ierr = VecGetArray(t, &tx);CHKERRQ(ierr);
 
1223
 
 
1224
    for (i = 0; i<n; i++) {
 
1225
      dpx += sx[i]*PetscConj(tx[i]);
 
1226
      nmx += tx[i]*PetscConj(tx[i]);
 
1227
    }
 
1228
    work[0] = dpx;
 
1229
    work[1] = nmx;
 
1230
    ierr = MPI_Allreduce(&work,&sum,2,MPIU_SCALAR,MPIU_SUM,((PetscObject)s)->comm);CHKERRQ(ierr);
 
1231
    *dp  = sum[0];
 
1232
    *nm  = sum[1];
 
1233
 
 
1234
    ierr = VecRestoreArray(t, &tx);CHKERRQ(ierr);
 
1235
    ierr = VecRestoreArray(s, &sx);CHKERRQ(ierr);
 
1236
    ierr = PetscLogFlops(4.0*n);CHKERRQ(ierr);
 
1237
  }
 
1238
  ierr = PetscLogEventBarrierEnd(VEC_DotNormBarrier,s,t,0,0,((PetscObject)s)->comm);CHKERRQ(ierr);
 
1239
  PetscFunctionReturn(0);
 
1240
}
 
1241
 
 
1242
#undef __FUNCT__  
 
1243
#define __FUNCT__ "VecSum"
 
1244
/*@
 
1245
   VecSum - Computes the sum of all the components of a vector.
 
1246
 
 
1247
   Collective on Vec
 
1248
 
 
1249
   Input Parameter:
 
1250
.  v - the vector 
 
1251
 
 
1252
   Output Parameter:
 
1253
.  sum - the result
 
1254
 
 
1255
   Level: beginner
 
1256
 
 
1257
   Concepts: sum^of vector entries
 
1258
 
 
1259
.seealso: VecNorm()
 
1260
@*/
 
1261
PetscErrorCode PETSCVEC_DLLEXPORT VecSum(Vec v,PetscScalar *sum)
 
1262
{
 
1263
  PetscErrorCode ierr;
 
1264
  PetscInt       i,n;
 
1265
  PetscScalar    *x,lsum = 0.0;
 
1266
 
 
1267
  PetscFunctionBegin;
 
1268
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
1269
  PetscValidScalarPointer(sum,2);
 
1270
  ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
1271
  ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
1272
  for (i=0; i<n; i++) {
 
1273
    lsum += x[i];
 
1274
  }
 
1275
  ierr = MPI_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,((PetscObject)v)->comm);CHKERRQ(ierr);
 
1276
  ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
1277
  PetscFunctionReturn(0);
 
1278
}
 
1279
 
 
1280
#undef __FUNCT__  
 
1281
#define __FUNCT__ "VecShift"
 
1282
/*@
 
1283
   VecShift - Shifts all of the components of a vector by computing
 
1284
   x[i] = x[i] + shift.
 
1285
 
 
1286
   Collective on Vec
 
1287
 
 
1288
   Input Parameters:
 
1289
+  v - the vector 
 
1290
-  shift - the shift
 
1291
 
 
1292
   Output Parameter:
 
1293
.  v - the shifted vector 
 
1294
 
 
1295
   Level: intermediate
 
1296
 
 
1297
   Concepts: vector^adding constant
 
1298
 
 
1299
@*/
 
1300
PetscErrorCode PETSCVEC_DLLEXPORT VecShift(Vec v,PetscScalar shift)
 
1301
{
 
1302
  PetscErrorCode ierr;
 
1303
  PetscInt       i,n;
 
1304
  PetscScalar    *x;
 
1305
 
 
1306
  PetscFunctionBegin;
 
1307
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
1308
  if (v->ops->shift) {
 
1309
    ierr = (*v->ops->shift)(v);CHKERRQ(ierr);
 
1310
  } else {
 
1311
    ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 
 
1312
    ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
1313
    for (i=0; i<n; i++) {
 
1314
      x[i] += shift;
 
1315
    }
 
1316
    ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
1317
  }
 
1318
  PetscFunctionReturn(0);
 
1319
}
 
1320
 
 
1321
#undef __FUNCT__  
 
1322
#define __FUNCT__ "VecAbs"
 
1323
/*@
 
1324
   VecAbs - Replaces every element in a vector with its absolute value.
 
1325
 
 
1326
   Collective on Vec
 
1327
 
 
1328
   Input Parameters:
 
1329
.  v - the vector 
 
1330
 
 
1331
   Level: intermediate
 
1332
 
 
1333
   Concepts: vector^absolute value
 
1334
 
 
1335
@*/
 
1336
PetscErrorCode PETSCVEC_DLLEXPORT VecAbs(Vec v)
 
1337
{
 
1338
  PetscErrorCode ierr;
 
1339
  PetscInt       i,n;
 
1340
  PetscScalar    *x;
 
1341
 
 
1342
  PetscFunctionBegin;
 
1343
  PetscValidHeaderSpecific(v,VEC_COOKIE,1);
 
1344
  if (v->ops->abs) {
 
1345
    ierr = (*v->ops->abs)(v);CHKERRQ(ierr);
 
1346
  } else {
 
1347
    ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr);
 
1348
    ierr = VecGetArray(v,&x);CHKERRQ(ierr);
 
1349
    for (i=0; i<n; i++) {
 
1350
      x[i] = PetscAbsScalar(x[i]);
 
1351
    }
 
1352
    ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
 
1353
  }
 
1354
  PetscFunctionReturn(0);
 
1355
}
 
1356
 
 
1357
#undef __FUNCT__
 
1358
#define __FUNCT__ "VecPermute"
 
1359
/*@
 
1360
  VecPermute - Permutes a vector in place using the given ordering.
 
1361
 
 
1362
  Input Parameters:
 
1363
+ vec   - The vector
 
1364
. order - The ordering
 
1365
- inv   - The flag for inverting the permutation
 
1366
 
 
1367
  Level: beginner
 
1368
 
 
1369
  Note: This function does not yet support parallel Index Sets
 
1370
 
 
1371
.seealso: MatPermute()
 
1372
.keywords: vec, permute
 
1373
@*/
 
1374
PetscErrorCode PETSCVEC_DLLEXPORT VecPermute(Vec x, IS row, PetscTruth inv)
 
1375
{
 
1376
  PetscScalar    *array, *newArray;
 
1377
  const PetscInt *idx;
 
1378
  PetscInt       i;
 
1379
  PetscErrorCode ierr;
 
1380
 
 
1381
  PetscFunctionBegin;
 
1382
  ierr = ISGetIndices(row, &idx);CHKERRQ(ierr);
 
1383
  ierr = VecGetArray(x, &array);CHKERRQ(ierr);
 
1384
  ierr = PetscMalloc(x->map->n*sizeof(PetscScalar), &newArray);CHKERRQ(ierr);
 
1385
#ifdef PETSC_USE_DEBUG
 
1386
  for(i = 0; i < x->map->n; i++) {
 
1387
    if ((idx[i] < 0) || (idx[i] >= x->map->n)) {
 
1388
      SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Permutation index %D is out of bounds: %D", i, idx[i]);
 
1389
    }
 
1390
  }
 
1391
#endif
 
1392
  if (!inv) {
 
1393
    for(i = 0; i < x->map->n; i++) newArray[i]      = array[idx[i]];
 
1394
  } else {
 
1395
    for(i = 0; i < x->map->n; i++) newArray[idx[i]] = array[i];
 
1396
  }
 
1397
  ierr = VecRestoreArray(x, &array);CHKERRQ(ierr);
 
1398
  ierr = ISRestoreIndices(row, &idx);CHKERRQ(ierr);
 
1399
  ierr = VecReplaceArray(x, newArray);CHKERRQ(ierr);
 
1400
  PetscFunctionReturn(0);
 
1401
}
 
1402
 
 
1403
#undef __FUNCT__  
 
1404
#define __FUNCT__ "VecEqual"
 
1405
/*@
 
1406
   VecEqual - Compares two vectors.
 
1407
 
 
1408
   Collective on Vec
 
1409
 
 
1410
   Input Parameters:
 
1411
+  vec1 - the first vector
 
1412
-  vec2 - the second vector
 
1413
 
 
1414
   Output Parameter:
 
1415
.  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.
 
1416
 
 
1417
   Level: intermediate
 
1418
 
 
1419
   Concepts: equal^two vectors
 
1420
   Concepts: vector^equality
 
1421
 
 
1422
@*/
 
1423
PetscErrorCode PETSCVEC_DLLEXPORT VecEqual(Vec vec1,Vec vec2,PetscTruth *flg)
 
1424
{
 
1425
  PetscScalar    *v1,*v2; 
 
1426
  PetscErrorCode ierr; 
 
1427
  PetscInt       n1,n2,N1,N2; 
 
1428
  PetscInt       state1,state2; 
 
1429
  PetscTruth     flg1; 
 
1430
 
 
1431
  PetscFunctionBegin; 
 
1432
  PetscValidHeaderSpecific(vec1,VEC_COOKIE,1); 
 
1433
  PetscValidHeaderSpecific(vec2,VEC_COOKIE,2); 
 
1434
  PetscValidPointer(flg,3); 
 
1435
  if (vec1 == vec2) { 
 
1436
    *flg = PETSC_TRUE; 
 
1437
  } else { 
 
1438
    ierr = VecGetSize(vec1,&N1);CHKERRQ(ierr); 
 
1439
    ierr = VecGetSize(vec2,&N2);CHKERRQ(ierr); 
 
1440
    if (N1 != N2) { 
 
1441
      flg1 = PETSC_FALSE; 
 
1442
    } else { 
 
1443
      ierr = VecGetLocalSize(vec1,&n1);CHKERRQ(ierr); 
 
1444
      ierr = VecGetLocalSize(vec2,&n2);CHKERRQ(ierr); 
 
1445
      if (n1 != n2) { 
 
1446
        flg1 = PETSC_FALSE; 
 
1447
      } else { 
 
1448
        ierr = PetscObjectStateQuery((PetscObject) vec1,&state1);CHKERRQ(ierr); 
 
1449
        ierr = PetscObjectStateQuery((PetscObject) vec2,&state2);CHKERRQ(ierr); 
 
1450
        ierr = VecGetArray(vec1,&v1);CHKERRQ(ierr); 
 
1451
        ierr = VecGetArray(vec2,&v2);CHKERRQ(ierr);
 
1452
#if defined(PETSC_USE_COMPLEX)
 
1453
        PetscInt k;
 
1454
        for (k=0; k<n1; k++){
 
1455
          if (PetscRealPart(v1[k]) != PetscRealPart(v2[k]) || PetscImaginaryPart(v1[k]) != PetscImaginaryPart(v2[k])){
 
1456
            flg1 = PETSC_FALSE;
 
1457
            break;
 
1458
          }
 
1459
        }
 
1460
#else 
 
1461
        ierr = PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);CHKERRQ(ierr); 
 
1462
#endif
 
1463
        ierr = VecRestoreArray(vec1,&v1);CHKERRQ(ierr); 
 
1464
        ierr = VecRestoreArray(vec2,&v2);CHKERRQ(ierr); 
 
1465
        ierr = PetscObjectSetState((PetscObject) vec1,state1);CHKERRQ(ierr); 
 
1466
        ierr = PetscObjectSetState((PetscObject) vec2,state2);CHKERRQ(ierr); 
 
1467
      } 
 
1468
    } 
 
1469
    /* combine results from all processors */ 
 
1470
    ierr = MPI_Allreduce(&flg1,flg,1,MPI_INT,MPI_MIN,((PetscObject)vec1)->comm);CHKERRQ(ierr); 
 
1471
  } 
 
1472
  PetscFunctionReturn(0); 
 
1473
}
 
1474
 
 
1475
 
 
1476