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 <stdio.h>
36 #include <stdlib.h>
37 #ifdef HAVE_MALLOC_H
38 #include <malloc.h>
39 #endif
40 #include <string.h>
41 #include <math.h>
42 #ifdef _OPENMP
43 #include <omp.h>
44 #endif
45 #ifdef USE_MPI
46 #include <mpi.h>
47 #endif
48 #include "lislib.h"
49
50 /************************************************
51 * lis_matrix_set_blocksize
52 * lis_matrix_convert
53 * lis_matrix_convert_self
54 * lis_matrix_copy
55 * lis_matrix_copyDLU
56 * lis_matrix_axpy
57 * lis_matrix_xpay
58 * lis_matrix_axpyz
59 * lis_matrix_scale
60 * lis_matrix_get_diagonal
61 * lis_matrix_shift_diagonal
62 * lis_matrix_shift_matrix
63 * lis_matrix_split
64 * lis_matrix_merge
65 * lis_matrix_solve
66 * lis_matrix_solveh
67 ************************************************/
68
69 #undef __FUNC__
70 #define __FUNC__ "lis_matrix_set_blocksize"
lis_matrix_set_blocksize(LIS_MATRIX A,LIS_INT bnr,LIS_INT bnc,LIS_INT row[],LIS_INT col[])71 LIS_INT lis_matrix_set_blocksize(LIS_MATRIX A, LIS_INT bnr, LIS_INT bnc, LIS_INT row[], LIS_INT col[])
72 {
73 LIS_INT i,n;
74 LIS_INT err;
75 LIS_INT *conv_row,*conv_col;
76
77 LIS_DEBUG_FUNC_IN;
78
79 err = lis_matrix_check(A,LIS_MATRIX_CHECK_NULL);
80 if( err ) return err;
81
82 if( bnr<=0 || bnc<=0 )
83 {
84 LIS_SETERR2(LIS_ERR_ILL_ARG,"bnr=%D <= 0 or bnc=%D <= 0\n",bnr,bnc);
85 return LIS_ERR_ILL_ARG;
86 }
87 if( (row==NULL && col!=NULL) || (row!=NULL && col==NULL) )
88 {
89 LIS_SETERR2(LIS_ERR_ILL_ARG,"either row[]=%x or col[]=%x is NULL\n",row,col);
90 return LIS_ERR_ILL_ARG;
91 }
92 if( row==NULL )
93 {
94 A->conv_bnr = bnr;
95 A->conv_bnc = bnc;
96 }
97 else
98 {
99 n = A->n;
100 conv_row = (LIS_INT *)lis_malloc(n*sizeof(LIS_INT),"lis_matrix_set_blocksize::conv_row");
101 if( conv_row==NULL )
102 {
103 LIS_SETERR_MEM(n*sizeof(LIS_INT));
104 return LIS_OUT_OF_MEMORY;
105 }
106 conv_col = (LIS_INT *)lis_malloc(n*sizeof(LIS_INT),"lis_matrix_set_blocksize::conv_col");
107 if( conv_col==NULL )
108 {
109 LIS_SETERR_MEM(n*sizeof(LIS_INT));
110 lis_free(conv_row);
111 return LIS_OUT_OF_MEMORY;
112 }
113 for(i=0;i<n;i++)
114 {
115 conv_row[i] = row[i];
116 conv_col[i] = col[i];
117 }
118 A->conv_row = conv_row;
119 A->conv_col = conv_col;
120 }
121
122 LIS_DEBUG_FUNC_OUT;
123 return LIS_SUCCESS;
124 }
125
126 #undef __FUNC__
127 #define __FUNC__ "lis_matrix_convert"
lis_matrix_convert(LIS_MATRIX Ain,LIS_MATRIX Aout)128 LIS_INT lis_matrix_convert(LIS_MATRIX Ain, LIS_MATRIX Aout)
129 {
130 LIS_INT err;
131 LIS_INT istmp;
132 LIS_INT convert_matrix_type;
133 LIS_MATRIX Atmp,Atmp2;
134
135 LIS_DEBUG_FUNC_IN;
136
137 err = lis_matrix_check(Ain,LIS_MATRIX_CHECK_ALL);
138 if( err ) return err;
139 err = lis_matrix_check(Aout,LIS_MATRIX_CHECK_NULL);
140 if( err ) return err;
141
142 err = lis_matrix_merge(Ain);
143 if( err ) return err;
144
145 convert_matrix_type = Aout->matrix_type;
146
147 if( Ain->matrix_type==convert_matrix_type && !Ain->is_block )
148 {
149 err = lis_matrix_copy(Ain,Aout);
150 return err;
151 }
152 if( Ain->matrix_type!=LIS_MATRIX_CSR )
153 {
154 istmp = LIS_TRUE;
155 switch( Ain->matrix_type )
156 {
157 case LIS_MATRIX_RCO:
158 switch( convert_matrix_type )
159 {
160 case LIS_MATRIX_CSR:
161 err = lis_matrix_convert_rco2csr(Ain,Aout);
162 LIS_DEBUG_FUNC_OUT;
163 return err;
164 break;
165 case LIS_MATRIX_BSR:
166 err = lis_matrix_convert_rco2bsr(Ain,Aout);
167 LIS_DEBUG_FUNC_OUT;
168 return err;
169 case LIS_MATRIX_CSC:
170 err = lis_matrix_convert_rco2csc(Ain,Aout);
171 LIS_DEBUG_FUNC_OUT;
172 return err;
173 break;
174 default:
175 err = lis_matrix_duplicate(Ain,&Atmp);
176 if( err ) return err;
177 err = lis_matrix_convert_rco2csr(Ain,Atmp);
178 break;
179 }
180 break;
181 case LIS_MATRIX_CSC:
182 switch( convert_matrix_type )
183 {
184 case LIS_MATRIX_BSC:
185 err = lis_matrix_convert_csc2bsc(Ain,Aout);
186 LIS_DEBUG_FUNC_OUT;
187 return err;
188 default:
189 err = lis_matrix_duplicate(Ain,&Atmp);
190 if( err ) return err;
191 err = lis_matrix_convert_csc2csr(Ain,Atmp);
192 break;
193 }
194 break;
195 case LIS_MATRIX_MSR:
196 err = lis_matrix_duplicate(Ain,&Atmp);
197 if( err ) return err;
198 err = lis_matrix_convert_msr2csr(Ain,Atmp);
199 break;
200 case LIS_MATRIX_DIA:
201 err = lis_matrix_duplicate(Ain,&Atmp);
202 if( err ) return err;
203 err = lis_matrix_convert_dia2csr(Ain,Atmp);
204 break;
205 case LIS_MATRIX_ELL:
206 err = lis_matrix_duplicate(Ain,&Atmp);
207 if( err ) return err;
208 err = lis_matrix_convert_ell2csr(Ain,Atmp);
209 break;
210 case LIS_MATRIX_JAD:
211 err = lis_matrix_duplicate(Ain,&Atmp);
212 if( err ) return err;
213 err = lis_matrix_convert_jad2csr(Ain,Atmp);
214 break;
215 case LIS_MATRIX_BSR:
216 err = lis_matrix_duplicate(Ain,&Atmp);
217 if( err ) return err;
218 err = lis_matrix_convert_bsr2csr(Ain,Atmp);
219 break;
220 case LIS_MATRIX_BSC:
221 err = lis_matrix_duplicate(Ain,&Atmp);
222 if( err ) return err;
223 err = lis_matrix_convert_bsc2csr(Ain,Atmp);
224 break;
225 case LIS_MATRIX_VBR:
226 err = lis_matrix_duplicate(Ain,&Atmp);
227 if( err ) return err;
228 err = lis_matrix_convert_vbr2csr(Ain,Atmp);
229 break;
230 case LIS_MATRIX_DNS:
231 err = lis_matrix_duplicate(Ain,&Atmp);
232 if( err ) return err;
233 err = lis_matrix_convert_dns2csr(Ain,Atmp);
234 break;
235 case LIS_MATRIX_COO:
236 err = lis_matrix_duplicate(Ain,&Atmp);
237 if( err ) return err;
238 err = lis_matrix_convert_coo2csr(Ain,Atmp);
239 break;
240 default:
241 LIS_SETERR_IMP;
242 err = LIS_ERR_NOT_IMPLEMENTED;
243 }
244 if( err )
245 {
246 return err;
247 }
248 if( convert_matrix_type== LIS_MATRIX_CSR )
249 {
250 lis_matrix_storage_destroy(Aout);
251 lis_matrix_DLU_destroy(Aout);
252 lis_matrix_diag_destroy(Aout->WD);
253 if( Aout->l2g_map ) lis_free( Aout->l2g_map );
254 if( Aout->commtable ) lis_commtable_destroy( Aout->commtable );
255 if( Aout->ranges ) lis_free( Aout->ranges );
256 lis_matrix_copy_struct(Atmp,Aout);
257 lis_free(Atmp);
258 return LIS_SUCCESS;
259 }
260 }
261 else
262 {
263 istmp = LIS_FALSE;
264 Atmp = Ain;
265 }
266
267 switch( convert_matrix_type )
268 {
269 case LIS_MATRIX_BSR:
270 err = lis_matrix_convert_csr2bsr(Atmp,Aout);
271 break;
272 case LIS_MATRIX_CSC:
273 err = lis_matrix_convert_csr2csc(Atmp,Aout);
274 break;
275 case LIS_MATRIX_MSR:
276 err = lis_matrix_convert_csr2msr(Atmp,Aout);
277 break;
278 case LIS_MATRIX_ELL:
279 err = lis_matrix_convert_csr2ell(Atmp,Aout);
280 break;
281 case LIS_MATRIX_DIA:
282 err = lis_matrix_convert_csr2dia(Atmp,Aout);
283 break;
284 case LIS_MATRIX_JAD:
285 err = lis_matrix_convert_csr2jad(Atmp,Aout);
286 break;
287 case LIS_MATRIX_BSC:
288 err = lis_matrix_duplicate(Atmp,&Atmp2);
289 if( err ) return err;
290 err = lis_matrix_convert_csr2csc(Atmp,Atmp2);
291 if( err ) return err;
292 if( Atmp!=Ain )
293 {
294 lis_matrix_destroy(Atmp);
295 }
296 Atmp = Atmp2;
297 istmp = LIS_TRUE;
298 err = lis_matrix_convert_csc2bsc(Atmp,Aout);
299 break;
300 case LIS_MATRIX_VBR:
301 err = lis_matrix_convert_csr2vbr(Atmp,Aout);
302 break;
303 case LIS_MATRIX_DNS:
304 err = lis_matrix_convert_csr2dns(Atmp,Aout);
305 break;
306 case LIS_MATRIX_COO:
307 err = lis_matrix_convert_csr2coo(Atmp,Aout);
308 break;
309 default:
310 LIS_SETERR_IMP;
311 err = LIS_ERR_NOT_IMPLEMENTED;
312 }
313
314 if( istmp ) lis_matrix_destroy(Atmp);
315 if( err )
316 {
317 return err;
318 }
319
320 LIS_DEBUG_FUNC_OUT;
321 return err;
322 }
323
324 #undef __FUNC__
325 #define __FUNC__ "lis_matrix_convert_self"
lis_matrix_convert_self(LIS_SOLVER solver)326 LIS_INT lis_matrix_convert_self(LIS_SOLVER solver)
327 {
328 LIS_INT err;
329 LIS_INT storage,block;
330 LIS_MATRIX A,B;
331
332 LIS_DEBUG_FUNC_IN;
333
334 A = solver->A;
335 storage = solver->options[LIS_OPTIONS_STORAGE];
336 block = solver->options[LIS_OPTIONS_STORAGE_BLOCK];
337
338 if( storage>0 && A->matrix_type!=storage )
339 {
340 err = lis_matrix_duplicate(A,&B);
341 if( err ) return err;
342 lis_matrix_set_blocksize(B,block,block,NULL,NULL);
343 lis_matrix_set_type(B,storage);
344 err = lis_matrix_convert(A,B);
345 if( err ) return err;
346 lis_matrix_storage_destroy(A);
347 lis_matrix_DLU_destroy(A);
348 lis_matrix_diag_destroy(A->WD);
349 if( A->l2g_map ) lis_free( A->l2g_map );
350 if( A->commtable ) lis_commtable_destroy( A->commtable );
351 if( A->ranges ) lis_free( A->ranges );
352 err = lis_matrix_copy_struct(B,A);
353 if( err ) return err;
354 lis_free(B);
355 if( A->matrix_type==LIS_MATRIX_JAD )
356 {
357 A->work = (LIS_SCALAR *)lis_malloc(A->n*sizeof(LIS_SCALAR),"lis_precon_create_bjacobi::A->work");
358 if( A->work==NULL )
359 {
360 LIS_SETERR_MEM(A->n*sizeof(LIS_SCALAR));
361 return LIS_OUT_OF_MEMORY;
362 }
363 }
364 }
365
366 LIS_DEBUG_FUNC_OUT;
367 return LIS_SUCCESS;
368 }
369
370 #undef __FUNC__
371 #define __FUNC__ "lis_matrix_copy"
lis_matrix_copy(LIS_MATRIX Ain,LIS_MATRIX Aout)372 LIS_INT lis_matrix_copy(LIS_MATRIX Ain, LIS_MATRIX Aout)
373 {
374 LIS_INT err;
375
376 LIS_DEBUG_FUNC_IN;
377
378 err = lis_matrix_check(Ain,LIS_MATRIX_CHECK_ALL);
379 if( err ) return err;
380 err = lis_matrix_check(Ain,LIS_MATRIX_CHECK_NULL);
381 if( err ) return err;
382
383 switch( Ain->matrix_type )
384 {
385 case LIS_MATRIX_CSR:
386 err = lis_matrix_copy_csr(Ain,Aout);
387 break;
388 case LIS_MATRIX_CSC:
389 err = lis_matrix_copy_csc(Ain,Aout);
390 break;
391 case LIS_MATRIX_MSR:
392 err = lis_matrix_copy_msr(Ain,Aout);
393 break;
394 case LIS_MATRIX_DIA:
395 err = lis_matrix_copy_dia(Ain,Aout);
396 break;
397 case LIS_MATRIX_ELL:
398 err = lis_matrix_copy_ell(Ain,Aout);
399 break;
400 case LIS_MATRIX_JAD:
401 err = lis_matrix_copy_jad(Ain,Aout);
402 break;
403 case LIS_MATRIX_BSR:
404 err = lis_matrix_copy_bsr(Ain,Aout);
405 break;
406 case LIS_MATRIX_VBR:
407 err = lis_matrix_copy_vbr(Ain,Aout);
408 break;
409 case LIS_MATRIX_DNS:
410 err = lis_matrix_copy_dns(Ain,Aout);
411 break;
412 case LIS_MATRIX_COO:
413 err = lis_matrix_copy_coo(Ain,Aout);
414 break;
415 case LIS_MATRIX_BSC:
416 err = lis_matrix_copy_bsc(Ain,Aout);
417 break;
418 default:
419 LIS_SETERR_IMP;
420 return LIS_ERR_NOT_IMPLEMENTED;
421 break;
422 }
423 LIS_DEBUG_FUNC_OUT;
424 return LIS_SUCCESS;
425 }
426
427 #undef __FUNC__
428 #define __FUNC__ "lis_matrix_copyDLU"
lis_matrix_copyDLU(LIS_MATRIX Ain,LIS_MATRIX_DIAG * D,LIS_MATRIX * L,LIS_MATRIX * U)429 LIS_INT lis_matrix_copyDLU(LIS_MATRIX Ain, LIS_MATRIX_DIAG *D, LIS_MATRIX *L, LIS_MATRIX *U)
430 {
431 LIS_INT err;
432
433 LIS_DEBUG_FUNC_IN;
434
435 err = lis_matrix_check(Ain,LIS_MATRIX_CHECK_ALL);
436 if( err ) return err;
437
438 switch( Ain->matrix_type )
439 {
440 case LIS_MATRIX_CSR:
441 err = lis_matrix_copyDLU_csr(Ain,D,L,U);
442 break;
443 /*
444 case LIS_MATRIX_CSC:
445 err = lis_matrix_copy_csc(Ain,Aout);
446 break;
447 case LIS_MATRIX_MSR:
448 err = lis_matrix_copy_msr(Ain,Aout);
449 break;
450 case LIS_MATRIX_DIA:
451 err = lis_matrix_copy_dia(Ain,Aout);
452 break;
453 case LIS_MATRIX_ELL:
454 err = lis_matrix_copy_ell(Ain,Aout);
455 break;
456 case LIS_MATRIX_JAD:
457 err = lis_matrix_copy_jad(Ain,Aout);
458 break;
459 case LIS_MATRIX_BSR:
460 err = lis_matrix_copy_bsr(Ain,Aout);
461 break;
462 case LIS_MATRIX_BSC:
463 err = lis_matrix_copy_bsc(Ain,Aout);
464 break;
465 case LIS_MATRIX_VBR:
466 err = lis_matrix_copy_vbr(Ain,Aout);
467 break;
468 case LIS_MATRIX_DNS:
469 err = lis_matrix_copy_dns(Ain,Aout);
470 break;
471 case LIS_MATRIX_COO:
472 err = lis_matrix_copy_coo(Ain,Aout);
473 break;
474 */
475 default:
476 LIS_SETERR_IMP;
477 *D = NULL;
478 *L = NULL;
479 *U = NULL;
480 return LIS_ERR_NOT_IMPLEMENTED;
481 break;
482 }
483 LIS_DEBUG_FUNC_OUT;
484 return LIS_SUCCESS;
485 }
486
487 #undef __FUNC__
488 #define __FUNC__ "lis_matrix_axpy"
lis_matrix_axpy(LIS_SCALAR alpha,LIS_MATRIX A,LIS_MATRIX B)489 LIS_INT lis_matrix_axpy(LIS_SCALAR alpha, LIS_MATRIX A, LIS_MATRIX B)
490 {
491 LIS_INT err;
492
493 LIS_DEBUG_FUNC_IN;
494
495 err = lis_matrix_check(A,LIS_MATRIX_CHECK_ALL);
496 if( err ) return err;
497
498 err = lis_matrix_check(B,LIS_MATRIX_CHECK_ALL);
499 if( err ) return err;
500
501 switch( A->matrix_type )
502 {
503 case LIS_MATRIX_DNS:
504 err = lis_matrix_axpy_dns(alpha,A,B);
505 break;
506
507 default:
508 LIS_SETERR_IMP;
509 return LIS_ERR_NOT_IMPLEMENTED;
510 break;
511 }
512 LIS_DEBUG_FUNC_OUT;
513 return LIS_SUCCESS;
514 }
515
516 #undef __FUNC__
517 #define __FUNC__ "lis_matrix_xpay"
lis_matrix_xpay(LIS_SCALAR alpha,LIS_MATRIX A,LIS_MATRIX B)518 LIS_INT lis_matrix_xpay(LIS_SCALAR alpha, LIS_MATRIX A, LIS_MATRIX B)
519 {
520 LIS_INT err;
521
522 LIS_DEBUG_FUNC_IN;
523
524 err = lis_matrix_check(A,LIS_MATRIX_CHECK_ALL);
525 if( err ) return err;
526
527 err = lis_matrix_check(B,LIS_MATRIX_CHECK_ALL);
528 if( err ) return err;
529
530 switch( A->matrix_type )
531 {
532 case LIS_MATRIX_DNS:
533 err = lis_matrix_xpay_dns(alpha,A,B);
534 break;
535
536 default:
537 LIS_SETERR_IMP;
538 return LIS_ERR_NOT_IMPLEMENTED;
539 break;
540 }
541 LIS_DEBUG_FUNC_OUT;
542 return LIS_SUCCESS;
543 }
544
545 #undef __FUNC__
546 #define __FUNC__ "lis_matrix_axpyz"
lis_matrix_axpyz(LIS_SCALAR alpha,LIS_MATRIX A,LIS_MATRIX B,LIS_MATRIX C)547 LIS_INT lis_matrix_axpyz(LIS_SCALAR alpha, LIS_MATRIX A, LIS_MATRIX B, LIS_MATRIX C)
548 {
549 LIS_INT err;
550
551 LIS_DEBUG_FUNC_IN;
552
553 err = lis_matrix_check(A,LIS_MATRIX_CHECK_ALL);
554 if( err ) return err;
555
556 err = lis_matrix_check(B,LIS_MATRIX_CHECK_ALL);
557 if( err ) return err;
558
559 err = lis_matrix_check(C,LIS_MATRIX_CHECK_ALL);
560 if( err ) return err;
561
562 switch( A->matrix_type )
563 {
564 case LIS_MATRIX_DNS:
565 err = lis_matrix_axpyz_dns(alpha,A,B,C);
566 break;
567
568 default:
569 LIS_SETERR_IMP;
570 return LIS_ERR_NOT_IMPLEMENTED;
571 break;
572 }
573 LIS_DEBUG_FUNC_OUT;
574 return LIS_SUCCESS;
575 }
576
577 #undef __FUNC__
578 #define __FUNC__ "lis_matrix_scale"
lis_matrix_scale(LIS_MATRIX A,LIS_VECTOR B,LIS_VECTOR D,LIS_INT action)579 LIS_INT lis_matrix_scale(LIS_MATRIX A, LIS_VECTOR B, LIS_VECTOR D, LIS_INT action)
580 {
581 LIS_INT i,n,np;
582 LIS_SCALAR *b,*d;
583
584 n = A->n;
585 np = A->np;
586 b = B->value;
587 d = D->value;
588
589 lis_matrix_get_diagonal(A,D);
590 if( action==LIS_SCALE_SYMM_DIAG )
591 {
592 #ifdef USE_MPI
593 if( A->np>D->np )
594 {
595 D->value = (LIS_SCALAR *)lis_realloc(D->value,A->np*sizeof(LIS_SCALAR));
596 if( D->value==NULL )
597 {
598 LIS_SETERR_MEM(A->np*sizeof(LIS_SCALAR));
599 return LIS_OUT_OF_MEMORY;
600 }
601 d = D->value;
602 }
603 lis_send_recv(A->commtable,d);
604 #endif
605 #ifdef _OPENMP
606 #pragma omp parallel for private(i)
607 #endif
608 for(i=0; i<np; i++)
609 {
610 d[i] = 1.0 / sqrt(fabs(d[i]));
611 }
612
613 switch( A->matrix_type )
614 {
615 case LIS_MATRIX_CSR:
616 lis_matrix_scale_symm_csr(A, d);
617 break;
618 case LIS_MATRIX_CSC:
619 lis_matrix_scale_symm_csc(A, d);
620 break;
621 case LIS_MATRIX_MSR:
622 lis_matrix_scale_symm_msr(A, d);
623 break;
624 case LIS_MATRIX_DIA:
625 lis_matrix_scale_symm_dia(A, d);
626 break;
627 case LIS_MATRIX_ELL:
628 lis_matrix_scale_symm_ell(A, d);
629 break;
630 case LIS_MATRIX_JAD:
631 lis_matrix_scale_symm_jad(A, d);
632 break;
633 case LIS_MATRIX_BSR:
634 lis_matrix_scale_symm_bsr(A, d);
635 break;
636 case LIS_MATRIX_BSC:
637 lis_matrix_scale_symm_bsc(A, d);
638 break;
639 case LIS_MATRIX_DNS:
640 lis_matrix_scale_symm_dns(A, d);
641 break;
642 case LIS_MATRIX_COO:
643 lis_matrix_scale_symm_coo(A, d);
644 break;
645 case LIS_MATRIX_VBR:
646 lis_matrix_scale_symm_vbr(A, d);
647 break;
648 default:
649 LIS_SETERR_IMP;
650 return LIS_ERR_NOT_IMPLEMENTED;
651 break;
652 }
653 }
654 else if (action==LIS_SCALE_JACOBI)
655 {
656 #ifdef _OPENMP
657 #pragma omp parallel for private(i)
658 #endif
659 for(i=0; i<n; i++)
660 {
661 d[i] = 1.0 / d[i];
662 }
663 switch( A->matrix_type )
664 {
665 case LIS_MATRIX_CSR:
666 lis_matrix_scale_csr(A, d);
667 break;
668 case LIS_MATRIX_CSC:
669 lis_matrix_scale_csc(A, d);
670 break;
671 case LIS_MATRIX_MSR:
672 lis_matrix_scale_msr(A, d);
673 break;
674 case LIS_MATRIX_DIA:
675 lis_matrix_scale_dia(A, d);
676 break;
677 case LIS_MATRIX_ELL:
678 lis_matrix_scale_ell(A, d);
679 break;
680 case LIS_MATRIX_JAD:
681 lis_matrix_scale_jad(A, d);
682 break;
683 case LIS_MATRIX_BSR:
684 lis_matrix_scale_bsr(A, d);
685 break;
686 case LIS_MATRIX_BSC:
687 lis_matrix_scale_bsc(A, d);
688 break;
689 case LIS_MATRIX_DNS:
690 lis_matrix_scale_dns(A, d);
691 break;
692 case LIS_MATRIX_COO:
693 lis_matrix_scale_coo(A, d);
694 break;
695 case LIS_MATRIX_VBR:
696 lis_matrix_scale_vbr(A, d);
697 break;
698 default:
699 LIS_SETERR_IMP;
700 return LIS_ERR_NOT_IMPLEMENTED;
701 break;
702 }
703 }
704
705 #ifdef _OPENMP
706 #pragma omp parallel for private(i)
707 #endif
708 for(i=0; i<n; i++)
709 {
710 b[i] = b[i]*d[i];
711 }
712 A->is_scaled = LIS_TRUE;
713 B->is_scaled = LIS_TRUE;
714 return LIS_SUCCESS;
715 }
716
717 /*NEH support for extended "solve_kernel" workflow*/
718 #undef __FUNC__
719 #define __FUNC__ "lis_matrix_psd_reset_scale"
lis_matrix_psd_reset_scale(LIS_MATRIX A)720 LIS_INT lis_matrix_psd_reset_scale(LIS_MATRIX A)
721 {
722 A->is_scaled=LIS_FALSE;
723 return LIS_SUCCESS;
724 }
725
726 #undef __FUNC__
727 #define __FUNC__ "lis_matrix_get_diagonal"
lis_matrix_get_diagonal(LIS_MATRIX A,LIS_VECTOR D)728 LIS_INT lis_matrix_get_diagonal(LIS_MATRIX A, LIS_VECTOR D)
729 {
730 LIS_SCALAR *d;
731
732 LIS_DEBUG_FUNC_IN;
733
734 d = D->value;
735 switch( A->matrix_type )
736 {
737 case LIS_MATRIX_CSR:
738 lis_matrix_get_diagonal_csr(A, d);
739 break;
740 case LIS_MATRIX_CSC:
741 lis_matrix_get_diagonal_csc(A, d);
742 break;
743 case LIS_MATRIX_MSR:
744 lis_matrix_get_diagonal_msr(A, d);
745 break;
746 case LIS_MATRIX_DIA:
747 lis_matrix_get_diagonal_dia(A, d);
748 break;
749 case LIS_MATRIX_ELL:
750 lis_matrix_get_diagonal_ell(A, d);
751 break;
752 case LIS_MATRIX_JAD:
753 lis_matrix_get_diagonal_jad(A, d);
754 break;
755 case LIS_MATRIX_BSR:
756 lis_matrix_get_diagonal_bsr(A, d);
757 break;
758 case LIS_MATRIX_BSC:
759 lis_matrix_get_diagonal_bsc(A, d);
760 break;
761 case LIS_MATRIX_DNS:
762 lis_matrix_get_diagonal_dns(A, d);
763 break;
764 case LIS_MATRIX_COO:
765 lis_matrix_get_diagonal_coo(A, d);
766 break;
767 case LIS_MATRIX_VBR:
768 lis_matrix_get_diagonal_vbr(A, d);
769 break;
770 default:
771 LIS_SETERR_IMP;
772 return LIS_ERR_NOT_IMPLEMENTED;
773 break;
774 }
775 LIS_DEBUG_FUNC_OUT;
776 return LIS_SUCCESS;
777 }
778
779 #undef __FUNC__
780 #define __FUNC__ "lis_matrix_shift_diagonal"
lis_matrix_shift_diagonal(LIS_MATRIX A,LIS_SCALAR sigma)781 LIS_INT lis_matrix_shift_diagonal(LIS_MATRIX A, LIS_SCALAR sigma)
782 {
783
784 LIS_DEBUG_FUNC_IN;
785
786 switch( A->matrix_type )
787 {
788 case LIS_MATRIX_CSR:
789 lis_matrix_shift_diagonal_csr(A, sigma);
790 break;
791 case LIS_MATRIX_CSC:
792 lis_matrix_shift_diagonal_csc(A, sigma);
793 break;
794 case LIS_MATRIX_MSR:
795 lis_matrix_shift_diagonal_msr(A, sigma);
796 break;
797 case LIS_MATRIX_DIA:
798 lis_matrix_shift_diagonal_dia(A, sigma);
799 break;
800 case LIS_MATRIX_ELL:
801 lis_matrix_shift_diagonal_ell(A, sigma);
802 break;
803 case LIS_MATRIX_JAD:
804 lis_matrix_shift_diagonal_jad(A, sigma);
805 break;
806 case LIS_MATRIX_BSR:
807 lis_matrix_shift_diagonal_bsr(A, sigma);
808 break;
809 case LIS_MATRIX_BSC:
810 lis_matrix_shift_diagonal_bsc(A, sigma);
811 break;
812 case LIS_MATRIX_DNS:
813 lis_matrix_shift_diagonal_dns(A, sigma);
814 break;
815 case LIS_MATRIX_COO:
816 lis_matrix_shift_diagonal_coo(A, sigma);
817 break;
818 case LIS_MATRIX_VBR:
819 lis_matrix_shift_diagonal_vbr(A, sigma);
820 break;
821
822 default:
823 LIS_SETERR_IMP;
824 return LIS_ERR_NOT_IMPLEMENTED;
825 break;
826 }
827 LIS_DEBUG_FUNC_OUT;
828 return LIS_SUCCESS;
829 }
830
831 #undef __FUNC__
832 #define __FUNC__ "lis_matrix_shift_matrix"
lis_matrix_shift_matrix(LIS_MATRIX A,LIS_MATRIX B,LIS_SCALAR sigma)833 LIS_INT lis_matrix_shift_matrix(LIS_MATRIX A, LIS_MATRIX B, LIS_SCALAR sigma)
834 {
835 LIS_INT err;
836 LIS_MATRIX Atmp,Btmp;
837
838 LIS_DEBUG_FUNC_IN;
839
840 err = lis_matrix_duplicate(A,&Atmp);
841 err = lis_matrix_duplicate(B,&Btmp);
842 if( err ) return err;
843 lis_matrix_set_type(Atmp,LIS_MATRIX_DNS);
844 err = lis_matrix_convert(A,Atmp);
845 if( err ) return err;
846 lis_matrix_set_type(Btmp,LIS_MATRIX_DNS);
847 err = lis_matrix_convert(B,Btmp);
848 if( err ) return err;
849 lis_matrix_axpy_dns(-sigma,Btmp,Atmp);
850 lis_matrix_convert(Atmp,A);
851 lis_matrix_destroy(Atmp);
852 lis_matrix_destroy(Btmp);
853
854 LIS_DEBUG_FUNC_OUT;
855 return LIS_SUCCESS;
856 }
857
858 #undef __FUNC__
859 #define __FUNC__ "lis_matrix_split"
lis_matrix_split(LIS_MATRIX A)860 LIS_INT lis_matrix_split(LIS_MATRIX A)
861 {
862 LIS_INT err;
863
864 LIS_DEBUG_FUNC_IN;
865
866 if( A->is_splited )
867 {
868 LIS_DEBUG_FUNC_OUT;
869 return LIS_SUCCESS;
870 }
871 switch( A->matrix_type )
872 {
873 case LIS_MATRIX_CSR:
874 err = lis_matrix_split_csr(A);
875 break;
876 case LIS_MATRIX_CSC:
877 err = lis_matrix_split_csc(A);
878 break;
879 case LIS_MATRIX_BSR:
880 err = lis_matrix_split_bsr(A);
881 break;
882 case LIS_MATRIX_MSR:
883 err = lis_matrix_split_msr(A);
884 break;
885 case LIS_MATRIX_ELL:
886 err = lis_matrix_split_ell(A);
887 break;
888 case LIS_MATRIX_DIA:
889 err = lis_matrix_split_dia(A);
890 break;
891 case LIS_MATRIX_JAD:
892 err = lis_matrix_split_jad(A);
893 break;
894 case LIS_MATRIX_BSC:
895 err = lis_matrix_split_bsc(A);
896 break;
897 case LIS_MATRIX_DNS:
898 err = lis_matrix_split_dns(A);
899 break;
900 case LIS_MATRIX_COO:
901 err = lis_matrix_split_coo(A);
902 break;
903 case LIS_MATRIX_VBR:
904 err = lis_matrix_split_vbr(A);
905 break;
906 default:
907 LIS_SETERR_IMP;
908 return LIS_ERR_NOT_IMPLEMENTED;
909 break;
910 }
911
912 if( err ) return err;
913 /*
914 if( !A->is_save )
915 {
916 lis_matrix_storage_destroy(A);
917 }
918 */
919 A->is_splited = LIS_TRUE;
920
921 LIS_DEBUG_FUNC_OUT;
922 return LIS_SUCCESS;
923 }
924
925 /*NEH support for extended "solve_kernel" workflow*/
926 #undef __FUNC__
927 #define __FUNC__ "lis_matrix_split_create"
lis_matrix_split_create(LIS_MATRIX A)928 LIS_INT lis_matrix_split_create(LIS_MATRIX A)
929 {
930 LIS_INT err;
931
932 LIS_DEBUG_FUNC_IN;
933
934 /* don't try to split if it is already split */
935 if( A->is_splited )
936 {
937 LIS_DEBUG_FUNC_OUT;
938 return LIS_SUCCESS;
939 }
940
941 switch( A->matrix_type )
942 {
943 case LIS_MATRIX_CSR:
944 err = lis_matrix_split_create_csr(A);
945 break;
946 case LIS_MATRIX_CSC:
947 LIS_SETERR_IMP;
948 return LIS_ERR_NOT_IMPLEMENTED;
949 case LIS_MATRIX_BSR:
950 LIS_SETERR_IMP;
951 return LIS_ERR_NOT_IMPLEMENTED;
952 case LIS_MATRIX_MSR:
953 LIS_SETERR_IMP;
954 return LIS_ERR_NOT_IMPLEMENTED;
955 case LIS_MATRIX_ELL:
956 LIS_SETERR_IMP;
957 return LIS_ERR_NOT_IMPLEMENTED;
958 case LIS_MATRIX_DIA:
959 LIS_SETERR_IMP;
960 return LIS_ERR_NOT_IMPLEMENTED;
961 case LIS_MATRIX_JAD:
962 LIS_SETERR_IMP;
963 return LIS_ERR_NOT_IMPLEMENTED;
964 case LIS_MATRIX_BSC:
965 LIS_SETERR_IMP;
966 return LIS_ERR_NOT_IMPLEMENTED;
967 case LIS_MATRIX_DNS:
968 LIS_SETERR_IMP;
969 return LIS_ERR_NOT_IMPLEMENTED;
970 case LIS_MATRIX_COO:
971 LIS_SETERR_IMP;
972 return LIS_ERR_NOT_IMPLEMENTED;
973 case LIS_MATRIX_VBR:
974 LIS_SETERR_IMP;
975 return LIS_ERR_NOT_IMPLEMENTED;
976 default:
977 LIS_SETERR_IMP;
978 return LIS_ERR_NOT_IMPLEMENTED;
979 }
980
981 if( err ) return err;
982
983 A->is_splited = LIS_TRUE;
984
985 LIS_DEBUG_FUNC_OUT;
986 return LIS_SUCCESS;
987 }
988
989 /*NEH support for extended "solve_kernel" workflow*/
990 #undef __FUNC__
991 #define __FUNC__ "lis_matrix_split_update"
lis_matrix_split_update(LIS_MATRIX A)992 LIS_INT lis_matrix_split_update(LIS_MATRIX A)
993 {
994 LIS_INT err;
995
996 LIS_DEBUG_FUNC_IN;
997
998 /* don't try to perform "split_update" if the matrix is not already split */
999 if( !A->is_splited )
1000 {
1001 return LIS_ERR_ILL_ARG;
1002 }
1003
1004 switch( A->matrix_type )
1005 {
1006 case LIS_MATRIX_CSR:
1007 err = lis_matrix_split_update_csr(A);
1008 break;
1009 case LIS_MATRIX_CSC:
1010 LIS_SETERR_IMP;
1011 return LIS_ERR_NOT_IMPLEMENTED;
1012 case LIS_MATRIX_BSR:
1013 LIS_SETERR_IMP;
1014 return LIS_ERR_NOT_IMPLEMENTED;
1015 case LIS_MATRIX_MSR:
1016 LIS_SETERR_IMP;
1017 return LIS_ERR_NOT_IMPLEMENTED;
1018 case LIS_MATRIX_ELL:
1019 LIS_SETERR_IMP;
1020 return LIS_ERR_NOT_IMPLEMENTED;
1021 case LIS_MATRIX_DIA:
1022 LIS_SETERR_IMP;
1023 return LIS_ERR_NOT_IMPLEMENTED;
1024 case LIS_MATRIX_JAD:
1025 LIS_SETERR_IMP;
1026 return LIS_ERR_NOT_IMPLEMENTED;
1027 case LIS_MATRIX_BSC:
1028 LIS_SETERR_IMP;
1029 return LIS_ERR_NOT_IMPLEMENTED;
1030 case LIS_MATRIX_DNS:
1031 LIS_SETERR_IMP;
1032 return LIS_ERR_NOT_IMPLEMENTED;
1033 case LIS_MATRIX_COO:
1034 LIS_SETERR_IMP;
1035 return LIS_ERR_NOT_IMPLEMENTED;
1036 case LIS_MATRIX_VBR:
1037 LIS_SETERR_IMP;
1038 return LIS_ERR_NOT_IMPLEMENTED;
1039 default:
1040 LIS_SETERR_IMP;
1041 return LIS_ERR_NOT_IMPLEMENTED;
1042 }
1043
1044 if( err ) return err;
1045
1046 LIS_DEBUG_FUNC_OUT;
1047 return LIS_SUCCESS;
1048 }
1049
1050 #undef __FUNC__
1051 #define __FUNC__ "lis_matrix_merge"
lis_matrix_merge(LIS_MATRIX A)1052 LIS_INT lis_matrix_merge(LIS_MATRIX A)
1053 {
1054 LIS_INT err;
1055
1056 LIS_DEBUG_FUNC_IN;
1057
1058 if( !A->is_splited || (A->is_save && A->is_splited) )
1059 {
1060 LIS_DEBUG_FUNC_OUT;
1061 return LIS_SUCCESS;
1062 }
1063 switch( A->matrix_type )
1064 {
1065 case LIS_MATRIX_CSR:
1066 err = lis_matrix_merge_csr(A);
1067 break;
1068 case LIS_MATRIX_CSC:
1069 err = lis_matrix_merge_csc(A);
1070 break;
1071 case LIS_MATRIX_MSR:
1072 err = lis_matrix_merge_msr(A);
1073 break;
1074 case LIS_MATRIX_BSR:
1075 err = lis_matrix_merge_bsr(A);
1076 break;
1077 case LIS_MATRIX_ELL:
1078 err = lis_matrix_merge_ell(A);
1079 break;
1080 case LIS_MATRIX_JAD:
1081 err = lis_matrix_merge_jad(A);
1082 break;
1083 case LIS_MATRIX_DIA:
1084 err = lis_matrix_merge_dia(A);
1085 break;
1086 case LIS_MATRIX_BSC:
1087 err = lis_matrix_merge_bsc(A);
1088 break;
1089 case LIS_MATRIX_DNS:
1090 err = lis_matrix_merge_dns(A);
1091 break;
1092 case LIS_MATRIX_COO:
1093 err = lis_matrix_merge_coo(A);
1094 break;
1095 case LIS_MATRIX_VBR:
1096 err = lis_matrix_merge_vbr(A);
1097 break;
1098 default:
1099 LIS_SETERR_IMP;
1100 return LIS_ERR_NOT_IMPLEMENTED;
1101 break;
1102 }
1103
1104 if( err ) return err;
1105 if( !A->is_save )
1106 {
1107 lis_matrix_DLU_destroy(A);
1108 A->is_splited = LIS_FALSE;
1109 }
1110
1111
1112 LIS_DEBUG_FUNC_OUT;
1113 return LIS_SUCCESS;
1114 }
1115
1116 #undef __FUNC__
1117 #define __FUNC__ "lis_matrix_solve"
lis_matrix_solve(LIS_MATRIX A,LIS_VECTOR b,LIS_VECTOR x,LIS_INT flag)1118 LIS_INT lis_matrix_solve(LIS_MATRIX A, LIS_VECTOR b, LIS_VECTOR x, LIS_INT flag)
1119 {
1120 LIS_DEBUG_FUNC_IN;
1121
1122 if( !A->is_splited ) lis_matrix_split(A);
1123
1124 switch( A->matrix_type )
1125 {
1126 case LIS_MATRIX_CSR:
1127 lis_matrix_solve_csr(A,b,x,flag);
1128 break;
1129 case LIS_MATRIX_BSR:
1130 lis_matrix_solve_bsr(A,b,x,flag);
1131 break;
1132 case LIS_MATRIX_CSC:
1133 lis_matrix_solve_csc(A,b,x,flag);
1134 break;
1135 case LIS_MATRIX_MSR:
1136 lis_matrix_solve_msr(A,b,x,flag);
1137 break;
1138 case LIS_MATRIX_ELL:
1139 lis_matrix_solve_ell(A,b,x,flag);
1140 break;
1141 case LIS_MATRIX_JAD:
1142 lis_matrix_solve_jad(A,b,x,flag);
1143 break;
1144 case LIS_MATRIX_DIA:
1145 lis_matrix_solve_dia(A,b,x,flag);
1146 break;
1147 case LIS_MATRIX_DNS:
1148 lis_matrix_solve_dns(A,b,x,flag);
1149 break;
1150 case LIS_MATRIX_BSC:
1151 lis_matrix_solve_bsc(A,b,x,flag);
1152 break;
1153 case LIS_MATRIX_VBR:
1154 lis_matrix_solve_vbr(A,b,x,flag);
1155 break;
1156 default:
1157 LIS_SETERR_IMP;
1158 return LIS_ERR_NOT_IMPLEMENTED;
1159 break;
1160 }
1161
1162 LIS_DEBUG_FUNC_OUT;
1163 return LIS_SUCCESS;
1164 }
1165
1166 #undef __FUNC__
1167 #define __FUNC__ "lis_matrix_solveh"
lis_matrix_solveh(LIS_MATRIX A,LIS_VECTOR b,LIS_VECTOR x,LIS_INT flag)1168 LIS_INT lis_matrix_solveh(LIS_MATRIX A, LIS_VECTOR b, LIS_VECTOR x, LIS_INT flag)
1169 {
1170 LIS_DEBUG_FUNC_IN;
1171
1172 if( !A->is_splited ) lis_matrix_split(A);
1173
1174 switch( A->matrix_type )
1175 {
1176 case LIS_MATRIX_CSR:
1177 lis_matrix_solveh_csr(A,b,x,flag);
1178 break;
1179 case LIS_MATRIX_BSR:
1180 lis_matrix_solveh_bsr(A,b,x,flag);
1181 break;
1182 case LIS_MATRIX_CSC:
1183 lis_matrix_solveh_csc(A,b,x,flag);
1184 break;
1185 case LIS_MATRIX_MSR:
1186 lis_matrix_solveh_msr(A,b,x,flag);
1187 break;
1188 case LIS_MATRIX_ELL:
1189 lis_matrix_solveh_ell(A,b,x,flag);
1190 break;
1191 case LIS_MATRIX_JAD:
1192 lis_matrix_solveh_jad(A,b,x,flag);
1193 break;
1194 case LIS_MATRIX_DIA:
1195 lis_matrix_solveh_dia(A,b,x,flag);
1196 break;
1197 case LIS_MATRIX_DNS:
1198 lis_matrix_solveh_dns(A,b,x,flag);
1199 break;
1200 case LIS_MATRIX_BSC:
1201 lis_matrix_solveh_bsc(A,b,x,flag);
1202 break;
1203 case LIS_MATRIX_VBR:
1204 lis_matrix_solveh_vbr(A,b,x,flag);
1205 break;
1206 default:
1207 LIS_SETERR_IMP;
1208 return LIS_ERR_NOT_IMPLEMENTED;
1209 break;
1210 }
1211
1212 LIS_DEBUG_FUNC_OUT;
1213 return LIS_SUCCESS;
1214 }
1215
1216
1217
1218
1219
1220