1 /* Copyright (C) 2005 The Scalable Software Infrastructure Project. All rights reserved.
2
3 Redistribution and use in source and binary forms, with or without
4 modification, are permitted provided that the following conditions are met:
5 1. Redistributions of source code must retain the above copyright
6 notice, this list of conditions and the following disclaimer.
7 2. Redistributions in binary form must reproduce the above copyright
8 notice, this list of conditions and the following disclaimer in the
9 documentation and/or other materials provided with the distribution.
10 3. Neither the name of the project nor the names of its contributors
11 may be used to endorse or promote products derived from this software
12 without specific prior written permission.
13
14 THIS SOFTWARE IS PROVIDED BY THE SCALABLE SOFTWARE INFRASTRUCTURE PROJECT
15 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
16 TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
17 PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE SCALABLE SOFTWARE INFRASTRUCTURE
18 PROJECT BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
19 OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
20 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
21 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
22 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
23 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
24 POSSIBILITY OF SUCH DAMAGE.
25 */
26
27 #ifdef HAVE_CONFIG_H
28 #include "lis_config.h"
29 #else
30 #ifdef HAVE_CONFIG_WIN_H
31 #include "lis_config_win.h"
32 #endif
33 #endif
34
35 #include <math.h>
36 #ifdef _OPENMP
37 #include <omp.h>
38 #endif
39 #ifdef USE_MPI
40 #include <mpi.h>
41 #endif
42 #include "lislib.h"
43
44 /********************************************************
45 * lis_vector_swap x <-> y
46 * lis_vector_copy y <- x
47 * lis_vector_axpy y <- y + alpha * x
48 * lis_vector_xpay y <- x + alpha * y
49 * lis_vector_axpyz z <- y + alpha * x
50 * lis_vector_scale x <- alpha * x
51 * lis_vector_pmul z_i <- x_i * y_i
52 * lis_vector_pdiv z_i <- x_i / y_i
53 * lis_vector_set_all x_i <- alpha
54 * lis_vector_abs x_i <- |x_i|
55 * lis_vector_reciprocal x_i <- 1 / x_i
56 * lis_vector_conjugate x_i <- conj(x_i)
57 * lis_vector_shift x_i <- x_i - sigma
58 * lis_vector_cgs classical Gram-Schmidt
59 ********************************************************/
60
61
62 /********************/
63 /* x <-> y */
64 /********************/
65 #undef __FUNC__
66 #define __FUNC__ "lis_vector_swap"
lis_vector_swap(LIS_VECTOR vx,LIS_VECTOR vy)67 LIS_INT lis_vector_swap(LIS_VECTOR vx, LIS_VECTOR vy)
68 {
69 LIS_INT i,n;
70 LIS_SCALAR *x,*y,t;
71
72 LIS_DEBUG_FUNC_IN;
73
74 n = vx->n;
75 #ifndef NO_ERROR_CHECK
76 if( n!=vy->n )
77 {
78 LIS_SETERR(LIS_ERR_ILL_ARG,"length of vector x and y is not equal\n");
79 return LIS_ERR_ILL_ARG;
80 }
81 #endif
82
83 x = vx->value;
84 y = vy->value;
85 #ifdef USE_VEC_COMP
86 #pragma cdir nodep
87 #pragma _NEC ivdep
88 #endif
89 #ifdef _OPENMP
90 #pragma omp parallel for private(i)
91 #endif
92 for(i=0; i<n; i++)
93 {
94 t = y[i];
95 y[i] = x[i];
96 x[i] = t;
97 }
98
99 LIS_DEBUG_FUNC_OUT;
100 return LIS_SUCCESS;
101 }
102
103
104 /********************/
105 /* y <- x */
106 /********************/
107 #undef __FUNC__
108 #define __FUNC__ "lis_vector_copy"
lis_vector_copy(LIS_VECTOR vx,LIS_VECTOR vy)109 LIS_INT lis_vector_copy(LIS_VECTOR vx, LIS_VECTOR vy)
110 {
111 LIS_INT i,n;
112 LIS_SCALAR *x,*y;
113
114 LIS_DEBUG_FUNC_IN;
115
116 n = vx->n;
117 #ifndef NO_ERROR_CHECK
118 if( n!=vy->n )
119 {
120 LIS_SETERR(LIS_ERR_ILL_ARG,"length of vector x and y is not equal\n");
121 return LIS_ERR_ILL_ARG;
122 }
123 #endif
124
125 x = vx->value;
126 y = vy->value;
127 #ifdef USE_VEC_COMP
128 #pragma cdir nodep
129 #pragma _NEC ivdep
130 #endif
131 #ifdef _OPENMP
132 #pragma omp parallel for private(i)
133 #endif
134 for(i=0; i<n; i++)
135 {
136 y[i] = x[i];
137 }
138
139 LIS_DEBUG_FUNC_OUT;
140 return LIS_SUCCESS;
141 }
142
143
144 /**********************/
145 /* y <- y + alpha * x */
146 /**********************/
147 #undef __FUNC__
148 #define __FUNC__ "lis_vector_axpy"
lis_vector_axpy(LIS_SCALAR alpha,LIS_VECTOR vx,LIS_VECTOR vy)149 LIS_INT lis_vector_axpy(LIS_SCALAR alpha, LIS_VECTOR vx, LIS_VECTOR vy)
150 {
151 LIS_INT i,n;
152 LIS_SCALAR *x,*y;
153
154 LIS_DEBUG_FUNC_IN;
155
156 n = vx->n;
157 #ifndef NO_ERROR_CHECK
158 if( n!=vy->n )
159 {
160 LIS_SETERR(LIS_ERR_ILL_ARG,"length of vector x and y is not equal\n");
161 return LIS_ERR_ILL_ARG;
162 }
163 #endif
164
165 x = vx->value;
166 y = vy->value;
167 #ifdef USE_VEC_COMP
168 #pragma cdir nodep
169 #pragma _NEC ivdep
170 #endif
171 #ifdef _OPENMP
172 #pragma omp parallel for private(i)
173 #endif
174 for(i=0; i<n; i++)
175 {
176 y[i] += alpha * x[i];
177 }
178
179 LIS_DEBUG_FUNC_OUT;
180 return LIS_SUCCESS;
181 }
182
183
184 /**********************/
185 /* y <- x + alpha * y */
186 /**********************/
187 #undef __FUNC__
188 #define __FUNC__ "lis_vector_xpay"
lis_vector_xpay(LIS_VECTOR vx,LIS_SCALAR alpha,LIS_VECTOR vy)189 LIS_INT lis_vector_xpay(LIS_VECTOR vx, LIS_SCALAR alpha, LIS_VECTOR vy)
190 {
191 LIS_INT i,n;
192 LIS_SCALAR *x,*y;
193
194 LIS_DEBUG_FUNC_IN;
195
196 n = vx->n;
197 #ifndef NO_ERROR_CHECK
198 if( n!=vy->n )
199 {
200 LIS_SETERR(LIS_ERR_ILL_ARG,"length of vector x and y is not equal\n");
201 return LIS_ERR_ILL_ARG;
202 }
203 #endif
204
205 x = vx->value;
206 y = vy->value;
207 #ifdef _OPENMP
208 #pragma omp parallel for private(i)
209 #endif
210 #ifdef USE_VEC_COMP
211 #pragma cdir nodep
212 #pragma _NEC ivdep
213 #endif
214 for(i=0; i<n; i++)
215 {
216 y[i] = x[i] + alpha * y[i];
217 }
218 LIS_DEBUG_FUNC_OUT;
219 return LIS_SUCCESS;
220 }
221
222
223 /**********************/
224 /* z <- y + alpha * x */
225 /**********************/
226 #undef __FUNC__
227 #define __FUNC__ "lis_vector_axpyz"
lis_vector_axpyz(LIS_SCALAR alpha,LIS_VECTOR vx,LIS_VECTOR vy,LIS_VECTOR vz)228 LIS_INT lis_vector_axpyz(LIS_SCALAR alpha, LIS_VECTOR vx, LIS_VECTOR vy, LIS_VECTOR vz)
229 {
230 LIS_INT i,n;
231 LIS_SCALAR *x,*y,*z;
232
233 LIS_DEBUG_FUNC_IN;
234
235 n = vx->n;
236 #ifndef NO_ERROR_CHECK
237 if( n!=vy->n || n!=vz->n )
238 {
239 LIS_SETERR(LIS_ERR_ILL_ARG,"length of vector x and y and z is not equal\n");
240 return LIS_ERR_ILL_ARG;
241 }
242 #endif
243
244 x = vx->value;
245 y = vy->value;
246 z = vz->value;
247 #ifdef USE_VEC_COMP
248 #pragma cdir nodep
249 #pragma _NEC ivdep
250 #endif
251 #ifdef _OPENMP
252 #pragma omp parallel for private(i)
253 #endif
254 for(i=0; i<n; i++)
255 {
256 z[i] = alpha * x[i] + y[i];
257 }
258 LIS_DEBUG_FUNC_OUT;
259 return LIS_SUCCESS;
260 }
261
262
263 /********************/
264 /* y <- alpha * x */
265 /********************/
266 #undef __FUNC__
267 #define __FUNC__ "lis_vector_scale"
lis_vector_scale(LIS_SCALAR alpha,LIS_VECTOR vx)268 LIS_INT lis_vector_scale(LIS_SCALAR alpha, LIS_VECTOR vx)
269 {
270 LIS_INT i,n;
271 LIS_SCALAR *x;
272
273 LIS_DEBUG_FUNC_IN;
274
275 n = vx->n;
276
277 x = vx->value;
278 #ifdef USE_VEC_COMP
279 #pragma cdir nodep
280 #pragma _NEC ivdep
281 #endif
282 #ifdef _OPENMP
283 #pragma omp parallel for private(i)
284 #endif
285 for(i=0; i<n; i++)
286 {
287 x[i] = alpha * x[i];
288 }
289
290 LIS_DEBUG_FUNC_OUT;
291 return LIS_SUCCESS;
292 }
293
294
295 /********************/
296 /* z_i <- x_i * y_i */
297 /********************/
298 #undef __FUNC__
299 #define __FUNC__ "lis_vector_pmul"
lis_vector_pmul(LIS_VECTOR vx,LIS_VECTOR vy,LIS_VECTOR vz)300 LIS_INT lis_vector_pmul(LIS_VECTOR vx,LIS_VECTOR vy,LIS_VECTOR vz)
301 {
302 LIS_INT i,n;
303 LIS_SCALAR *x,*y,*z;
304
305 LIS_DEBUG_FUNC_IN;
306
307 n = vx->n;
308 #ifndef NO_ERROR_CHECK
309 if( n!=vy->n || n!=vz->n )
310 {
311 LIS_SETERR(LIS_ERR_ILL_ARG,"length of vector x and y and z is not equal\n");
312 return LIS_ERR_ILL_ARG;
313 }
314 #endif
315
316 x = vx->value;
317 y = vy->value;
318 z = vz->value;
319 #ifdef USE_VEC_COMP
320 #pragma cdir nodep
321 #pragma _NEC ivdep
322 #endif
323 #ifdef _OPENMP
324 #pragma omp parallel for private(i)
325 #endif
326 for(i=0; i<n; i++)
327 {
328 z[i] = x[i] * y[i];
329 }
330 LIS_DEBUG_FUNC_OUT;
331 return LIS_SUCCESS;
332 }
333
334
335 /********************/
336 /* z_i <- x_i / y_i */
337 /********************/
338 #undef __FUNC__
339 #define __FUNC__ "lis_vector_pdiv"
lis_vector_pdiv(LIS_VECTOR vx,LIS_VECTOR vy,LIS_VECTOR vz)340 LIS_INT lis_vector_pdiv(LIS_VECTOR vx,LIS_VECTOR vy,LIS_VECTOR vz)
341 {
342 LIS_INT i,n;
343 LIS_SCALAR *x,*y,*z;
344
345 LIS_DEBUG_FUNC_IN;
346
347 n = vx->n;
348 #ifndef NO_ERROR_CHECK
349 if( n!=vy->n || n!=vz->n )
350 {
351 LIS_SETERR(LIS_ERR_ILL_ARG,"length of vector x and y and z is not equal\n");
352 return LIS_ERR_ILL_ARG;
353 }
354 #endif
355
356 x = vx->value;
357 y = vy->value;
358 z = vz->value;
359 #ifdef USE_VEC_COMP
360 #pragma cdir nodep
361 #pragma _NEC ivdep
362 #endif
363 #ifdef _OPENMP
364 #pragma omp parallel for private(i)
365 #endif
366 for(i=0; i<n; i++)
367 {
368 z[i] = x[i] / y[i];
369 }
370 LIS_DEBUG_FUNC_OUT;
371 return LIS_SUCCESS;
372 }
373
374
375 /********************/
376 /* x_i <- alpha */
377 /********************/
378 #undef __FUNC__
379 #define __FUNC__ "lis_vector_set_all"
lis_vector_set_all(LIS_SCALAR alpha,LIS_VECTOR vx)380 LIS_INT lis_vector_set_all(LIS_SCALAR alpha, LIS_VECTOR vx)
381 {
382 LIS_INT i,n;
383 LIS_SCALAR *x;
384
385 LIS_DEBUG_FUNC_IN;
386
387 n = vx->n;
388
389 x = vx->value;
390 #ifdef USE_VEC_COMP
391 #pragma cdir nodep
392 #pragma _NEC ivdep
393 #endif
394 #ifdef _OPENMP
395 #pragma omp parallel for private(i)
396 #endif
397 for(i=0; i<n; i++)
398 {
399 x[i] = alpha;
400 }
401
402 LIS_DEBUG_FUNC_OUT;
403 return LIS_SUCCESS;
404 }
405
406
407 /********************/
408 /* x_i <- |x_i| */
409 /********************/
410 #undef __FUNC__
411 #define __FUNC__ "lis_vector_abs"
lis_vector_abs(LIS_VECTOR vx)412 LIS_INT lis_vector_abs(LIS_VECTOR vx)
413 {
414 LIS_INT i,n;
415 LIS_SCALAR *x;
416
417 LIS_DEBUG_FUNC_IN;
418
419 x = vx->value;
420 n = vx->n;
421 #ifdef USE_VEC_COMP
422 #pragma cdir nodep
423 #pragma _NEC ivdep
424 #endif
425 #ifdef _OPENMP
426 #pragma omp parallel for private(i)
427 #endif
428 for(i=0; i<n; i++)
429 {
430 x[i] = fabs(x[i]);
431 }
432 LIS_DEBUG_FUNC_OUT;
433 return LIS_SUCCESS;
434 }
435
436
437 /********************/
438 /* x_i <- 1 / x_i */
439 /********************/
440 #undef __FUNC__
441 #define __FUNC__ "lis_vector_reciprocal"
lis_vector_reciprocal(LIS_VECTOR vx)442 LIS_INT lis_vector_reciprocal(LIS_VECTOR vx)
443 {
444 LIS_INT i,n;
445 LIS_SCALAR *x;
446
447 LIS_DEBUG_FUNC_IN;
448
449 x = vx->value;
450 n = vx->n;
451 #ifdef USE_VEC_COMP
452 #pragma cdir nodep
453 #pragma _NEC ivdep
454 #endif
455 #ifdef _OPENMP
456 #pragma omp parallel for private(i)
457 #endif
458 for(i=0; i<n; i++)
459 {
460 x[i] = 1.0 / x[i];
461 }
462 LIS_DEBUG_FUNC_OUT;
463 return LIS_SUCCESS;
464 }
465
466
467 /********************/
468 /* x_i <- conj(x_i) */
469 /********************/
470 #undef __FUNC__
471 #define __FUNC__ "lis_vector_conjugate"
lis_vector_conjugate(LIS_VECTOR vx)472 LIS_INT lis_vector_conjugate(LIS_VECTOR vx)
473 {
474 LIS_INT i,n;
475 LIS_SCALAR *x;
476
477 LIS_DEBUG_FUNC_IN;
478
479 x = vx->value;
480 n = vx->n;
481 #ifdef USE_VEC_COMP
482 #pragma cdir nodep
483 #pragma _NEC ivdep
484 #endif
485 #ifdef _OPENMP
486 #pragma omp parallel for private(i)
487 #endif
488 for(i=0; i<n; i++)
489 {
490 #ifdef _COMPLEX
491 x[i] = conj(x[i]);
492 #endif
493 }
494 LIS_DEBUG_FUNC_OUT;
495 return LIS_SUCCESS;
496 }
497
498
499 /************************/
500 /* x_i <- alpha + x_i */
501 /************************/
502 #undef __FUNC__
503 #define __FUNC__ "lis_vector_shift"
lis_vector_shift(LIS_SCALAR sigma,LIS_VECTOR vx)504 LIS_INT lis_vector_shift(LIS_SCALAR sigma, LIS_VECTOR vx)
505 {
506 LIS_INT i,n;
507 LIS_SCALAR *x;
508
509 LIS_DEBUG_FUNC_IN;
510
511 x = vx->value;
512 n = vx->n;
513 #ifdef USE_VEC_COMP
514 #pragma cdir nodep
515 #pragma _NEC ivdep
516 #endif
517 #ifdef _OPENMP
518 #pragma omp parallel for private(i)
519 #endif
520 for(i=0; i<n; i++)
521 {
522 x[i] = x[i] - sigma;
523 }
524 LIS_DEBUG_FUNC_OUT;
525 return LIS_SUCCESS;
526 }
527
528
529 /*************************************/
530 /* QR <- X by Classical Gram-Schmidt */
531 /*************************************/
532 #undef __FUNC__
533 #define __FUNC__ "lis_vector_cgs"
lis_vector_cgs(LIS_INT n,LIS_VECTOR * x,LIS_VECTOR * q,LIS_VECTOR * r)534 LIS_INT lis_vector_cgs(LIS_INT n, LIS_VECTOR *x, LIS_VECTOR *q, LIS_VECTOR *r)
535 {
536 LIS_INT i, j, k;
537 LIS_VECTOR x_k;
538 LIS_REAL nrm2;
539 LIS_REAL tol;
540
541 lis_vector_duplicate(x[0], &x_k);
542
543 tol = 1e-6;
544
545 for (k=0;k<n;k++)
546 {
547 lis_vector_set_all(0.0,q[k]);
548 lis_vector_set_all(0.0,r[k]);
549 }
550
551 for (k=0;k<n;k++)
552 {
553 lis_vector_copy(x[k],x_k);
554 for (j=0;j<k;j++)
555 {
556 r[k]->value[j] = 0;
557 for (i=0;i<n;i++)
558 {
559 r[k]->value[j] += q[j]->value[i] * x[k]->value[i];
560 }
561 for (i=0;i<n;i++)
562 {
563 x_k->value[i] += q[j]->value[i] * x[k]->value[i];
564 }
565 }
566 lis_vector_nrm2(x_k, &nrm2);
567 if (nrm2<tol) break;
568 for (i=0;i<n;i++)
569 {
570 q[k]->value[i] = x_k->value[i] / nrm2;
571 }
572 }
573
574 lis_vector_destroy(x_k);
575
576 return 0;
577 }
578
579