Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
blas_array.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 
8 #ifdef _OPENMP
9 #include <omp.h>
10 #endif
11 
12 #include "fasp.h"
13 #include "fasp_functs.h"
14 
15 /*---------------------------------*/
16 /*-- Public Functions --*/
17 /*---------------------------------*/
18 
35 void fasp_blas_array_ax (const INT n,
36  const REAL a,
37  REAL *x)
38 {
39  INT i;
40  INT use_openmp = FALSE;
41 
42 #ifdef _OPENMP
43  INT myid, mybegin, myend, nthreads;
44  if ( n > OPENMP_HOLDS ) {
45  use_openmp = TRUE;
46  nthreads = FASP_GET_NUM_THREADS();
47  }
48 #endif
49 
50  if (a == 1.0) {
51 
52  }
53  else {
54  if (use_openmp) {
55 #ifdef _OPENMP
56 #pragma omp parallel private(myid, mybegin, myend, i)
57  {
58  myid = omp_get_thread_num();
59  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
60  for (i=mybegin; i<myend; ++i) x[i] *= a;
61  }
62 #endif
63  }
64  else {
65  for (i=0; i<n; ++i) x[i] *= a;
66  }
67  }
68 }
69 
87 void fasp_blas_array_axpy (const INT n,
88  const REAL a,
89  REAL *x,
90  REAL *y)
91 {
92  INT i;
93  INT use_openmp = FALSE;
94 
95 #ifdef _OPENMP
96  INT myid, mybegin, myend, nthreads;
97  if ( n > OPENMP_HOLDS ) {
98  use_openmp = TRUE;
99  nthreads = FASP_GET_NUM_THREADS();
100  }
101 #endif
102 
103  if (a==1.0) {
104  if (use_openmp) {
105 #ifdef _OPENMP
106 #pragma omp parallel private(myid, mybegin, myend, i) num_threads(nthreads)
107  {
108  myid = omp_get_thread_num();
109  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
110  for (i=mybegin; i<myend; ++i) y[i] += x[i];
111  }
112 #endif
113  }
114  else {
115  for (i=0; i<n; ++i) y[i] += x[i];
116  }
117  }
118  else if (a==-1.0) {
119  if (use_openmp) {
120 #ifdef _OPENMP
121 #pragma omp parallel private(myid, mybegin, myend, i) num_threads(nthreads)
122  {
123  myid = omp_get_thread_num();
124  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
125  for (i=mybegin; i<myend; ++i) y[i] -= x[i];
126  }
127 #endif
128  }
129  else {
130  for (i=0; i<n; ++i) y[i] -= x[i];
131  }
132  }
133  else {
134  if (use_openmp) {
135 #ifdef _OPENMP
136 #pragma omp parallel private(myid, mybegin, myend, i) num_threads(nthreads)
137  {
138  myid = omp_get_thread_num();
139  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
140  for (i=mybegin; i<myend; ++i) y[i] += a*x[i];
141  }
142 #endif
143  }
144  else {
145  for (i=0; i<n; ++i) y[i] += a*x[i];
146  }
147  }
148 }
149 
168  const REAL a,
169  REAL *x,
170  REAL *y,
171  REAL *z)
172 {
173  INT i;
174  INT use_openmp = FALSE;
175 
176 #ifdef _OPENMP
177  INT myid, mybegin, myend, nthreads;
178  if ( n > OPENMP_HOLDS ) {
179  use_openmp = TRUE;
180  nthreads = FASP_GET_NUM_THREADS();
181  }
182 #endif
183 
184  if (use_openmp) {
185 #ifdef _OPENMP
186 #pragma omp parallel private(myid, mybegin, myend, i) num_threads(nthreads)
187  {
188  myid = omp_get_thread_num();
189  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
190  for (i=mybegin; i<myend; ++i) z[i] = a*x[i]+y[i];
191  }
192 #endif
193  }
194  else {
195  for (i=0; i<n; ++i) z[i] = a*x[i]+y[i];
196  }
197 }
198 
219  const REAL a,
220  REAL *x,
221  const REAL b,
222  REAL *y)
223 {
224  INT i;
225  INT use_openmp = FALSE;
226 
227 #ifdef _OPENMP
228  INT myid, mybegin, myend, nthreads;
229  if ( n > OPENMP_HOLDS ) {
230  use_openmp = TRUE;
231  nthreads = FASP_GET_NUM_THREADS();
232  }
233 #endif
234 
235  if (use_openmp) {
236 #ifdef _OPENMP
237 #pragma omp parallel private(myid, mybegin, myend, i) num_threads(nthreads)
238  {
239  myid = omp_get_thread_num();
240  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
241  for (i=mybegin; i<myend; ++i) y[i] = a*x[i]+b*y[i];
242  }
243 #endif
244  }
245  else {
246  for (i=0; i<n; ++i) y[i] = a*x[i]+b*y[i];
247  }
248 
249 }
250 
268  const REAL * x,
269  const REAL * y)
270 {
271  INT i;
272  INT use_openmp = FALSE;
273  REAL value = 0.0;
274 
275 #ifdef _OPENMP
276  if ( n > OPENMP_HOLDS ) use_openmp = TRUE;
277 #endif
278 
279  if (use_openmp) {
280 #ifdef _OPENMP
281 #pragma omp parallel for reduction(+:value) private(i)
282 #endif
283  for ( i=0; i<n; ++i ) value += x[i]*y[i];
284  }
285  else {
286  for ( i=0; i<n; ++i ) value += x[i]*y[i];
287  }
288 
289  return value;
290 }
291 
308  const REAL *x)
309 {
310  INT i;
311  INT use_openmp = FALSE;
312  REAL onenorm = 0.;
313 
314 #ifdef _OPENMP
315  if ( n > OPENMP_HOLDS ) {
316  use_openmp = TRUE;
317  }
318 #endif
319 
320  if (use_openmp) {
321 #ifdef _OPENMP
322 #pragma omp parallel for reduction(+:onenorm) private(i)
323 #endif
324  for (i=0;i<n;++i) onenorm+=ABS(x[i]);
325  }
326  else {
327  for (i=0;i<n;++i) onenorm+=ABS(x[i]);
328  }
329  return onenorm;
330 }
331 
348  const REAL *x)
349 {
350  INT i;
351  INT use_openmp = FALSE;
352  REAL twonorm = 0.;
353 
354 #ifdef _OPENMP
355  if ( n > OPENMP_HOLDS ) {
356  use_openmp = TRUE;
357  }
358 #endif
359 
360  if (use_openmp) {
361 #ifdef _OPENMP
362 #pragma omp parallel for reduction(+:twonorm) private(i)
363 #endif
364  for (i=0;i<n;++i) twonorm+=x[i]*x[i];
365  }
366  else {
367  for (i=0;i<n;++i) twonorm+=x[i]*x[i];
368  }
369 
370  return sqrt(twonorm);
371 }
372 
389  const REAL *x)
390 {
391  INT i;
392  INT use_openmp = FALSE;
393  REAL infnorm = 0.0;
394 
395 #ifdef _OPENMP
396  INT myid, mybegin, myend, nthreads;
397  if ( n > OPENMP_HOLDS ) {
398  use_openmp = TRUE;
399  nthreads = FASP_GET_NUM_THREADS();
400  }
401 #endif
402 
403  if (use_openmp) {
404 #ifdef _OPENMP
405  REAL infnorm_loc = 0.0;
406 #pragma omp parallel firstprivate(infnorm_loc) private(myid, mybegin, myend, i)
407  {
408  myid = omp_get_thread_num();
409  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
410  for (i=mybegin; i<myend; ++i) infnorm_loc = MAX(infnorm_loc, ABS(x[i]));
411  if (infnorm_loc > infnorm) {
412 #pragma omp critical
413  infnorm = MAX(infnorm_loc, infnorm);
414  }
415  }
416 #endif
417  }
418  else {
419  for (i=0;i<n;++i) infnorm=MAX(infnorm,ABS(x[i]));
420  }
421 
422  return infnorm;
423 }
424 
425 /*---------------------------------*/
426 /*-- End of File --*/
427 /*---------------------------------*/
#define TRUE
Definition of logic type.
Definition: fasp_const.h:67
#define REAL
Definition: fasp.h:67
REAL fasp_blas_array_norm1(const INT n, const REAL *x)
L1 norm of array x.
Definition: blas_array.c:307
#define INT
Definition: fasp.h:64
void fasp_blas_array_axpby(const INT n, const REAL a, REAL *x, const REAL b, REAL *y)
y = a*x + b*y
Definition: blas_array.c:218
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
#define OPENMP_HOLDS
Definition: fasp_const.h:248
Main header file for FASP.
void fasp_blas_array_axpy(const INT n, const REAL a, REAL *x, REAL *y)
y = a*x + y
Definition: blas_array.c:87
REAL fasp_blas_array_norminf(const INT n, const REAL *x)
Linf norm of array x.
Definition: blas_array.c:388
#define ABS(a)
Definition: fasp.h:74
REAL fasp_blas_array_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies (x,y)
Definition: blas_array.c:267
REAL fasp_blas_array_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: blas_array.c:347
void fasp_blas_array_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: blas_array.c:35
#define FALSE
Definition: fasp_const.h:68
void fasp_blas_array_axpyz(const INT n, const REAL a, REAL *x, REAL *y, REAL *z)
z = a*x + y
Definition: blas_array.c:167
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72