Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
blas_vec.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 
33 void fasp_blas_dvec_axpy (const REAL a,
34  dvector *x,
35  dvector *y)
36 {
37  INT i, m=x->row;
38  INT use_openmp = FALSE;
39  REAL *xpt=x->val, *ypt=y->val;
40 
41 #ifdef _OPENMP
42  INT myid, mybegin, myend, nthreads;
43  if ( m > OPENMP_HOLDS ) {
44  use_openmp = TRUE;
45  nthreads = FASP_GET_NUM_THREADS();
46  }
47 #endif
48 
49  if ((y->row-m)!=0) {
50  printf("### ERROR: Two vectors have different dimensions!\n");
51  fasp_chkerr(ERROR_DATA_STRUCTURE, __FUNCTION__);
52  }
53 
54  if (use_openmp) {
55 #ifdef _OPENMP
56 #pragma omp parallel private(myid,mybegin,myend,i) num_threads(nthreads)
57  {
58  myid = omp_get_thread_num();
59  FASP_GET_START_END(myid, nthreads, m, &mybegin, &myend);
60  for (i=mybegin; i<myend; ++i) ypt[i] += a*xpt[i];
61  }
62 #endif
63  }
64  else {
65  for (i=0; i<m; ++i) ypt[i] += a*xpt[i];
66  }
67 }
68 
85 void fasp_blas_dvec_axpyz(const REAL a,
86  dvector *x,
87  dvector *y,
88  dvector *z)
89 {
90  const INT m=x->row;
91  REAL *xpt=x->val, *ypt=y->val, *zpt=z->val;
92 
93  if ((y->row-m)!=0) {
94  printf("### ERROR: Two vectors have different dimensions!\n");
95  fasp_chkerr(ERROR_DATA_STRUCTURE, __FUNCTION__);
96  }
97 
98  z->row = m;
99 
100  memcpy(ypt, zpt, m*sizeof(dvector));
101  fasp_blas_array_axpy(m, a, xpt, zpt);
102 }
103 
122  dvector *y)
123 {
124  REAL value=0;
125  INT i;
126  const INT length=x->row;
127  REAL *xpt=x->val, *ypt=y->val;
128 
129  INT use_openmp = FALSE;
130 
131 #ifdef _OPENMP
132  if ( length > OPENMP_HOLDS ) {
133  use_openmp = TRUE;
134  }
135 #endif
136 
137  if (use_openmp) {
138 #ifdef _OPENMP
139 #pragma omp parallel for reduction(+:value) private(i)
140 #endif
141  for (i=0; i<length; ++i) value+=xpt[i]*ypt[i];
142  }
143  else {
144  for (i=0; i<length; ++i) value+=xpt[i]*ypt[i];
145  }
146  return value;
147 }
148 
149 
168  dvector *y)
169 {
170  REAL diff=0, temp=0;
171  INT i;
172  const INT length=x->row;
173  REAL *xpt=x->val, *ypt=y->val;
174 
175  INT use_openmp = FALSE;
176 
177  if (length!=y->row) {
178  printf("### ERROR: Two vectors have different dimensions!\n");
179  fasp_chkerr(ERROR_DATA_STRUCTURE, __FUNCTION__);
180  }
181 
182 #ifdef _OPENMP
183  if ( length > OPENMP_HOLDS ) {
184  use_openmp = TRUE;
185  }
186 #endif
187 
188  if (use_openmp) {
189 #ifdef _OPENMP
190 #pragma omp parallel for reduction(+:temp,diff) private(i)
191 #endif
192  for (i=0;i<length;++i) {
193  temp += xpt[i]*xpt[i];
194  diff += pow(xpt[i]-ypt[i],2);
195  }
196  }
197  else {
198  for (i=0;i<length;++i) {
199  temp += xpt[i]*xpt[i];
200  diff += pow(xpt[i]-ypt[i],2);
201  }
202  }
203  return sqrt(diff/temp);
204 }
205 
223 {
224  REAL onenorm=0;
225  INT i;
226  const INT length=x->row;
227  REAL *xpt=x->val;
228 
229  INT use_openmp = FALSE;
230 
231 #ifdef _OPENMP
232  if ( length > OPENMP_HOLDS ) {
233  use_openmp = TRUE;
234  }
235 #endif
236 
237  if (use_openmp) {
238 #ifdef _OPENMP
239 #pragma omp parallel for reduction(+:onenorm) private(i)
240 #endif
241  for (i=0;i<length;++i) onenorm+=ABS(xpt[i]);
242  }
243  else {
244  for (i=0;i<length;++i) onenorm+=ABS(xpt[i]);
245  }
246  return onenorm;
247 }
248 
266 {
267  REAL twonorm=0;
268  INT i;
269  const INT length=x->row;
270  REAL *xpt=x->val;
271 
272  INT use_openmp = FALSE;
273 
274 #ifdef _OPENMP
275  if ( length > OPENMP_HOLDS ) {
276  use_openmp = TRUE;
277  }
278 #endif
279 
280  if (use_openmp) {
281 #ifdef _OPENMP
282 #pragma omp parallel for reduction(+:twonorm) private(i)
283 #endif
284  for (i=0;i<length;++i) twonorm+=xpt[i]*xpt[i];
285  }
286  else {
287  for (i=0;i<length;++i) twonorm+=xpt[i]*xpt[i];
288  }
289  return sqrt(twonorm);
290 }
291 
306 {
307  INT i;
308  const INT length=x->row;
309  REAL *xpt=x->val;
310 
311  register REAL infnorm=0;
312 
313  for (i=0;i<length;++i) infnorm=MAX(infnorm,ABS(xpt[i]));
314 
315  return infnorm;
316 }
317 
318 /*---------------------------------*/
319 /*-- End of File --*/
320 /*---------------------------------*/
#define TRUE
Definition of logic type.
Definition: fasp_const.h:67
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
REAL * val
actual vector entries
Definition: fasp.h:348
Vector with n entries of REAL type.
Definition: fasp.h:342
#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
REAL fasp_blas_dvec_norm2(dvector *x)
L2 norm of dvector x.
Definition: blas_vec.c:265
REAL fasp_blas_dvec_norminf(dvector *x)
Linf norm of dvector x.
Definition: blas_vec.c:305
REAL fasp_blas_dvec_dotprod(dvector *x, dvector *y)
Inner product of two vectors (x,y)
Definition: blas_vec.c:121
#define OPENMP_HOLDS
Definition: fasp_const.h:248
REAL fasp_blas_dvec_relerr(dvector *x, dvector *y)
Relative error of two dvector x and y.
Definition: blas_vec.c:167
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
INT row
number of rows
Definition: fasp.h:345
#define ERROR_DATA_STRUCTURE
Definition: fasp_const.h:38
#define ABS(a)
Definition: fasp.h:74
void fasp_blas_dvec_axpy(const REAL a, dvector *x, dvector *y)
y = a*x + y
Definition: blas_vec.c:33
#define FALSE
Definition: fasp_const.h:68
void fasp_blas_dvec_axpyz(const REAL a, dvector *x, dvector *y, dvector *z)
z = a*x + y, z is a third vector (z is cleared)
Definition: blas_vec.c:85
REAL fasp_blas_dvec_norm1(dvector *x)
L1 norm of dvector x.
Definition: blas_vec.c:222
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72