Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
vec.c
Go to the documentation of this file.
1 
8 #include <math.h>
9 
10 #ifdef _OPENMP
11 #include <omp.h>
12 #endif
13 
14 #include "fasp.h"
15 #include "fasp_functs.h"
16 
17 /*---------------------------------*/
18 /*-- Public Functions --*/
19 /*---------------------------------*/
20 
34 {
35  INT i;
36 
37  for ( i=0; i<u->row; i++ ) {
38  if ( ISNAN(u->val[i]) ) return TRUE;
39  }
40 
41  return FALSE;
42 }
43 
57 {
58  dvector u;
59 
60  u.row = m;
61  u.val = (REAL *)fasp_mem_calloc(m,sizeof(REAL));
62 
63  return u;
64 }
65 
79 {
80  ivector u;
81 
82  u.row = m;
83  u.val = (INT *)fasp_mem_calloc(m,sizeof(INT));
84 
85  return u;
86 }
87 
99 void fasp_dvec_alloc (const INT m,
100  dvector *u)
101 {
102  u->row = m;
103  u->val = (REAL*)fasp_mem_calloc(m,sizeof(REAL));
104 
105  return;
106 }
107 
119 void fasp_ivec_alloc (const INT m,
120  ivector *u)
121 {
122 
123  u->row = m;
124  u->val = (INT*)fasp_mem_calloc(m,sizeof(INT));
125 
126  return;
127 }
128 
140 {
141  if (u==NULL) return;
142 
143  fasp_mem_free(u->val);
144  u->row = 0; u->val = NULL;
145 }
146 
160 {
161  if (u==NULL) return;
162 
163  fasp_mem_free(u->val);
164  u->row = 0; u->val = NULL;
165 }
166 
178 {
179  x->row = 0; x->val = NULL;
180 }
181 
203 void fasp_dvec_rand (const INT n,
204  dvector *x)
205 {
206  const INT va=(REAL) 0;
207  const INT vb=(REAL) n;
208 
209  INT s=1; srand(s);
210 
211  INT i,j;
212 
213  x->row = n;
214 
215  for (i=0; i<n; ++i){
216  j = 1 + (INT) (((REAL)n)*rand()/(RAND_MAX+1.0));
217  x->val[i] = (((REAL)j)-va)/(vb-va);
218  }
219 }
220 
236  dvector *x,
237  REAL val)
238 {
239  INT i;
240  REAL *xpt=x->val;
241 
242  if (n>0) x->row=n;
243  else n=x->row;
244 
245 #ifdef _OPENMP
246  // variables for OpenMP
247  INT myid, mybegin, myend;
248  INT nthreads = FASP_GET_NUM_THREADS();
249 #endif
250 
251  if (val == 0.0) {
252 
253 #ifdef _OPENMP
254  if (n > OPENMP_HOLDS) {
255 #pragma omp parallel for private(myid, mybegin, myend)
256  for (myid = 0; myid < nthreads; myid++ ) {
257  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
258  memset(&xpt[mybegin],0x0, sizeof(REAL)*(myend-mybegin));
259  }
260  }
261  else {
262 #endif
263  memset(xpt, 0x0, sizeof(REAL)*n);
264 #ifdef _OPENMP
265  }
266 #endif
267 
268  }
269 
270  else {
271 
272 #ifdef _OPENMP
273  if (n > OPENMP_HOLDS) {
274 #pragma omp parallel for private(myid, mybegin, myend)
275  for (myid = 0; myid < nthreads; myid++ ) {
276  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
277  for (i=mybegin; i<myend; ++i) xpt[i]=val;
278  }
279  }
280  else {
281 #endif
282  for (i=0; i<n; ++i) xpt[i]=val;
283 #ifdef _OPENMP
284  }
285 #endif
286 
287  }
288 }
289 
304 void fasp_ivec_set (const INT m,
305  ivector *u)
306 {
307  INT i;
308  INT n = u->row;
309  INT nthreads = 1, use_openmp = FALSE;
310 
311 #ifdef _OPENMP
312  if ( n > OPENMP_HOLDS ) {
313  use_openmp = TRUE;
314  nthreads = FASP_GET_NUM_THREADS();
315  }
316 #endif
317 
318  if (use_openmp) {
319  INT mybegin, myend, myid;
320 #ifdef _OPENMP
321 #pragma omp parallel for private(myid, mybegin, myend)
322 #endif
323  for (myid = 0; myid < nthreads; myid++ ) {
324  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
325  for (i=mybegin; i<myend; ++i) u->val[i] = m;
326  }
327  }
328  else {
329 
330  for (i=0; i<n; ++i) u->val[i] = m;
331  }
332 }
333 
346  dvector *y)
347 {
348  y->row=x->row;
349  memcpy(y->val,x->val,x->row*sizeof(REAL));
350 }
351 
369  dvector *y)
370 {
371  const INT length=x->row;
372  REAL Linf=0., diffi=0.;
373  REAL *xpt=x->val, *ypt=y->val;
374 
375  INT use_openmp = FALSE;
376  INT i;
377 
378 #ifdef _OPENMP
379  INT myid, mybegin, myend, nthreads;
380  if ( length > OPENMP_HOLDS ) {
381  use_openmp = TRUE;
382  nthreads = FASP_GET_NUM_THREADS();
383  }
384 #endif
385 
386  if(use_openmp) {
387 #ifdef _OPENMP
388  REAL temp = 0.;
389 #pragma omp parallel firstprivate(temp) private(myid, mybegin, myend, i, diffi)
390  {
391  myid = omp_get_thread_num();
392  FASP_GET_START_END(myid, nthreads, length, &mybegin, &myend);
393  for(i=mybegin; i<myend; i++) {
394  if ((diffi = ABS(xpt[i]-ypt[i])) > temp) temp = diffi;
395  }
396 #pragma omp critical
397  if(temp > Linf) Linf = temp;
398  }
399 #endif
400  }
401  else {
402  for (i=0; i<length; ++i) {
403  if ((diffi = ABS(xpt[i]-ypt[i])) > Linf) Linf = diffi;
404  }
405  }
406 
407  return Linf;
408 }
409 
422  dvector *diag)
423 {
424  // information about dvector
425  const INT n = b->row;
426  REAL *val = b->val;
427 
428  // local variables
429  INT i;
430 
431  if (diag->row != n) {
432  printf("### ERROR: Sizes of diag = %d != dvector = %d!", diag->row, n);
433  fasp_chkerr(ERROR_MISC, __FUNCTION__);
434  }
435 
436  INT use_openmp = FALSE;
437 
438 #ifdef _OPENMP
439  INT mybegin, myend, myid, nthreads;
440  if ( n > OPENMP_HOLDS ){
441  use_openmp = TRUE;
442  nthreads = FASP_GET_NUM_THREADS();
443  }
444 #endif
445 
446  if (use_openmp) {
447 #ifdef _OPENMP
448 #pragma omp parallel for private(myid, mybegin,myend)
449  for (myid = 0; myid < nthreads; myid++ ) {
450  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
451  for (i=mybegin; i<myend; ++i) val[i] = val[i]/sqrt(diag->val[i]);
452  }
453 #endif
454  }
455  else {
456  for (i=0; i<n; ++i) val[i] = val[i]/sqrt(diag->val[i]);
457  }
458 
459  return;
460 }
461 
462 /*---------------------------------*/
463 /*-- End of File --*/
464 /*---------------------------------*/
#define TRUE
Definition of logic type.
Definition: fasp_const.h:67
dvector fasp_dvec_create(const INT m)
Create dvector data space of REAL type.
Definition: vec.c:56
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: message.c:199
#define REAL
Definition: fasp.h:67
Vector with n entries of INT type.
Definition: fasp.h:356
void fasp_dvec_free(dvector *u)
Free vector data space of REAL type.
Definition: vec.c:139
void fasp_dvec_cp(dvector *x, dvector *y)
Copy dvector x to dvector y.
Definition: vec.c:345
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
Vector with n entries of REAL type.
Definition: fasp.h:342
INT fasp_dvec_isnan(dvector *u)
Check a dvector whether there is NAN.
Definition: vec.c:33
#define INT
Definition: fasp.h:64
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
#define ISNAN(a)
Definition: fasp.h:83
void fasp_ivec_free(ivector *u)
Free vector data space of INT type.
Definition: vec.c:159
void fasp_dvec_rand(const INT n, dvector *x)
Generate random REAL vector in the range from 0 to 1.
Definition: vec.c:203
#define OPENMP_HOLDS
Definition: fasp_const.h:248
#define ERROR_MISC
Definition: fasp_const.h:35
Main header file for FASP.
INT * val
actual vector entries
Definition: fasp.h:362
void fasp_ivec_set(const INT m, ivector *u)
Set ivector value to be m.
Definition: vec.c:304
INT row
number of rows
Definition: fasp.h:345
void fasp_dvec_symdiagscale(dvector *b, dvector *diag)
Symmetric diagonal scaling D^{-1/2}b.
Definition: vec.c:421
void fasp_ivec_alloc(const INT m, ivector *u)
Create vector data space of INT type.
Definition: vec.c:119
#define ABS(a)
Definition: fasp.h:74
REAL fasp_dvec_maxdiff(dvector *x, dvector *y)
Maximal difference of two dvector x and y.
Definition: vec.c:368
INT row
number of rows
Definition: fasp.h:359
void fasp_dvec_null(dvector *x)
Initialize dvector.
Definition: vec.c:177
void fasp_dvec_set(INT n, dvector *x, REAL val)
Initialize dvector x[i]=val for i=0:n-1.
Definition: vec.c:235
void fasp_dvec_alloc(const INT m, dvector *u)
Create dvector data space of REAL type.
Definition: vec.c:99
#define FALSE
Definition: fasp_const.h:68
ivector fasp_ivec_create(const INT m)
Create vector data space of INT type.
Definition: vec.c:78