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 <stdarg.h>
42 #include <math.h>
43 #ifdef _OPENMP
44 #include <omp.h>
45 #endif
46 #ifdef USE_MPI
47 #include <mpi.h>
48 #endif
49 #include "lislib.h"
50
51 /************************************************
52 * function | SOM |
53 *-----------------------------+-----+
54 * lis_matrix_set | o |
55 * lis_matrix_setDLU | o |
56 * lis_matrix_malloc | o |
57 * lis_matrix_elements_copy | o |
58 * lis_matrix_transpose | o |
59 * lis_matrix_split | o |
60 * lis_matrix_merge | o |
61 *-----------------------------+-----+-----+
62 * function |merge|split|
63 *-----------------------------+-----+-----|
64 * lis_matrix_convert | o | |
65 * lis_matrix_copy | o | o |
66 * lis_matrix_get_diagonal | o | o |
67 * lis_matrix_shift_diagonal | o | o |
68 * lis_matrix_scale | o | o |
69 * lis_matrix_scale_symm | o | o |
70 * lis_matrix_normf | o | o |
71 * lis_matrix_sort | o | o |
72 * lis_matrix_solve | xxx | o |
73 * lis_matrix_solveh | xxx | o |
74 ************************************************/
75
76 #undef __FUNC__
77 #define __FUNC__ "lis_matrix_set_coo"
lis_matrix_set_coo(LIS_INT nnz,LIS_INT * row,LIS_INT * col,LIS_SCALAR * value,LIS_MATRIX A)78 LIS_INT lis_matrix_set_coo(LIS_INT nnz, LIS_INT *row, LIS_INT *col, LIS_SCALAR *value, LIS_MATRIX A)
79 {
80 LIS_INT err;
81
82 LIS_DEBUG_FUNC_IN;
83
84 #if 0
85 err = lis_matrix_check(A,LIS_MATRIX_CHECK_SET);
86 if( err ) return err;
87 #else
88 if(lis_matrix_is_assembled(A)) {
89 LIS_DEBUG_FUNC_OUT;
90 return LIS_SUCCESS;
91 }
92 else {
93 err = lis_matrix_check(A,LIS_MATRIX_CHECK_SET);
94 if( err ) return err;
95 }
96 #endif
97
98 A->row = row;
99 A->col = col;
100 A->value = value;
101 A->is_copy = LIS_FALSE;
102 A->status = -LIS_MATRIX_COO;
103 A->nnz = nnz;
104
105 LIS_DEBUG_FUNC_OUT;
106 return LIS_SUCCESS;
107 }
108
109 #undef __FUNC__
110 #define __FUNC__ "lis_matrix_setDLU_coo"
lis_matrix_setDLU_coo(LIS_INT lnnz,LIS_INT unnz,LIS_SCALAR * diag,LIS_INT * lrow,LIS_INT * lcol,LIS_SCALAR * lvalue,LIS_INT * urow,LIS_INT * ucol,LIS_SCALAR * uvalue,LIS_MATRIX A)111 LIS_INT lis_matrix_setDLU_coo(LIS_INT lnnz, LIS_INT unnz, LIS_SCALAR *diag, LIS_INT *lrow, LIS_INT *lcol, LIS_SCALAR *lvalue, LIS_INT *urow, LIS_INT *ucol, LIS_SCALAR *uvalue, LIS_MATRIX A)
112 {
113 LIS_INT err;
114 LIS_MATRIX_DIAG D;
115
116 LIS_DEBUG_FUNC_IN;
117
118 #if 0
119 err = lis_matrix_check(A,LIS_MATRIX_CHECK_SET);
120 if( err ) return err;
121 #else
122 if(lis_matrix_is_assembled(A)) return LIS_SUCCESS;
123 else {
124 err = lis_matrix_check(A,LIS_MATRIX_CHECK_SET);
125 if( err ) return err;
126 }
127 #endif
128
129 A->L = (LIS_MATRIX_CORE)lis_calloc(sizeof(struct LIS_MATRIX_CORE_STRUCT),"lis_matrix_setDLU_coo::A->L");
130 if( A->L==NULL )
131 {
132 LIS_SETERR_MEM(sizeof(struct LIS_MATRIX_CORE_STRUCT));
133 return LIS_OUT_OF_MEMORY;
134 }
135 A->U = (LIS_MATRIX_CORE)lis_calloc(sizeof(struct LIS_MATRIX_CORE_STRUCT),"lis_matrix_setDLU_coo::A->U");
136 if( A->U==NULL )
137 {
138 LIS_SETERR_MEM(sizeof(struct LIS_MATRIX_CORE_STRUCT));
139 lis_matrix_DLU_destroy(A);
140 return LIS_OUT_OF_MEMORY;
141 }
142 err = lis_matrix_diag_create(A->n,0,A->comm,&D);
143 if( err )
144 {
145 lis_matrix_DLU_destroy(A);
146 return err;
147 }
148
149 lis_free(D->value);
150 D->value = diag;
151 A->D = D;
152 A->L->nnz = lnnz;
153 A->L->row = lrow;
154 A->L->col = lcol;
155 A->L->value = lvalue;
156 A->U->nnz = unnz;
157 A->U->row = urow;
158 A->U->col = ucol;
159 A->U->value = uvalue;
160 A->is_copy = LIS_FALSE;
161 A->status = -LIS_MATRIX_COO;
162 A->is_splited = LIS_TRUE;
163
164 LIS_DEBUG_FUNC_OUT;
165 return LIS_SUCCESS;
166 }
167
168 #undef __FUNC__
169 #define __FUNC__ "lis_matrix_malloc_coo"
lis_matrix_malloc_coo(LIS_INT nnz,LIS_INT ** row,LIS_INT ** col,LIS_SCALAR ** value)170 LIS_INT lis_matrix_malloc_coo(LIS_INT nnz, LIS_INT **row, LIS_INT **col, LIS_SCALAR **value)
171 {
172 LIS_DEBUG_FUNC_IN;
173
174 *row = NULL;
175 *col = NULL;
176 *value = NULL;
177
178 *row = (LIS_INT *)lis_malloc( nnz*sizeof(LIS_INT),"lis_matrix_malloc_coo::row" );
179 if( *row==NULL )
180 {
181 LIS_SETERR_MEM(nnz*sizeof(LIS_INT));
182 lis_free2(3,*row,*col,*value);
183 return LIS_OUT_OF_MEMORY;
184 }
185 *col = (LIS_INT *)lis_malloc( nnz*sizeof(LIS_INT),"lis_matrix_malloc_coo::col" );
186 if( *col==NULL )
187 {
188 LIS_SETERR_MEM(nnz*sizeof(LIS_INT));
189 lis_free2(3,*row,*col,*value);
190 return LIS_OUT_OF_MEMORY;
191 }
192 *value = (LIS_SCALAR *)lis_malloc( nnz*sizeof(LIS_SCALAR),"lis_matrix_malloc_coo::value" );
193 if( *value==NULL )
194 {
195 LIS_SETERR_MEM(nnz*sizeof(LIS_SCALAR));
196 lis_free2(3,*row,*col,*value);
197 return LIS_OUT_OF_MEMORY;
198 }
199 LIS_DEBUG_FUNC_OUT;
200 return LIS_SUCCESS;
201 }
202
203
204 #undef __FUNC__
205 #define __FUNC__ "lis_matrix_elements_copy_coo"
lis_matrix_elements_copy_coo(LIS_INT nnz,LIS_INT * row,LIS_INT * col,LIS_SCALAR * value,LIS_INT * o_row,LIS_INT * o_col,LIS_SCALAR * o_value)206 LIS_INT lis_matrix_elements_copy_coo(LIS_INT nnz, LIS_INT *row, LIS_INT *col, LIS_SCALAR *value, LIS_INT *o_row, LIS_INT *o_col, LIS_SCALAR *o_value)
207 {
208 LIS_INT i;
209
210 LIS_DEBUG_FUNC_IN;
211
212
213 #ifdef _OPENMP
214 #pragma omp parallel for private(i)
215 #endif
216 for(i=0;i<nnz;i++)
217 {
218 o_row[i] = row[i];
219 o_col[i] = col[i];
220 o_value[i] = value[i];
221 }
222
223 LIS_DEBUG_FUNC_OUT;
224 return LIS_SUCCESS;
225 }
226
227 #undef __FUNC__
228 #define __FUNC__ "lis_matrix_copy_coo"
lis_matrix_copy_coo(LIS_MATRIX Ain,LIS_MATRIX Aout)229 LIS_INT lis_matrix_copy_coo(LIS_MATRIX Ain, LIS_MATRIX Aout)
230 {
231 LIS_INT err;
232 LIS_INT i,n,nnz,lnnz,unnz;
233 LIS_INT *row,*col;
234 LIS_INT *lrow,*lcol;
235 LIS_INT *urow,*ucol;
236 LIS_SCALAR *value,*lvalue,*uvalue,*diag;
237
238 LIS_DEBUG_FUNC_IN;
239
240 n = Ain->n;
241
242 if( Ain->is_splited )
243 {
244 lnnz = Ain->L->nnz;
245 unnz = Ain->U->nnz;
246 lrow = NULL;
247 lcol = NULL;
248 lvalue = NULL;
249 urow = NULL;
250 ucol = NULL;
251 uvalue = NULL;
252 diag = NULL;
253
254 err = lis_matrix_malloc_coo(lnnz,&lrow,&lcol,&lvalue);
255 if( err )
256 {
257 return err;
258 }
259 err = lis_matrix_malloc_coo(unnz,&urow,&ucol,&uvalue);
260 if( err )
261 {
262 lis_free2(7,diag,urow,lcol,urow,lcol,uvalue,lvalue);
263 return err;
264 }
265 diag = (LIS_SCALAR *)lis_malloc(n*sizeof(LIS_SCALAR),"lis_matrix_copy_coo::diag");
266 if( diag==NULL )
267 {
268 lis_free2(7,diag,urow,lcol,urow,lcol,uvalue,lvalue);
269 return err;
270 }
271
272 #ifdef _OPENMP
273 #pragma omp parallel for private(i)
274 #endif
275 for(i=0;i<n;i++)
276 {
277 diag[i] = Ain->D->value[i];
278 }
279 lis_matrix_elements_copy_coo(lnnz,Ain->L->row,Ain->L->col,Ain->L->value,lrow,lcol,lvalue);
280 lis_matrix_elements_copy_coo(unnz,Ain->U->row,Ain->U->col,Ain->U->value,urow,ucol,uvalue);
281
282 err = lis_matrix_setDLU_coo(lnnz,unnz,diag,lrow,lcol,lvalue,urow,ucol,uvalue,Aout);
283 if( err )
284 {
285 lis_free2(7,diag,urow,lcol,urow,lcol,uvalue,lvalue);
286 return err;
287 }
288 }
289 if( !Ain->is_splited || (Ain->is_splited && Ain->is_save) )
290 {
291 row = NULL;
292 col = NULL;
293 value = NULL;
294 nnz = Ain->nnz;
295 err = lis_matrix_malloc_coo(nnz,&row,&col,&value);
296 if( err )
297 {
298 return err;
299 }
300
301 lis_matrix_elements_copy_coo(nnz,Ain->row,Ain->col,Ain->value,row,col,value);
302
303 err = lis_matrix_set_coo(nnz,row,col,value,Aout);
304 if( err )
305 {
306 lis_free2(3,row,col,value);
307 return err;
308 }
309 }
310 err = lis_matrix_assemble(Aout);
311 if( err )
312 {
313 lis_matrix_storage_destroy(Aout);
314 return err;
315 }
316 LIS_DEBUG_FUNC_OUT;
317 return LIS_SUCCESS;
318 }
319
320 #undef __FUNC__
321 #define __FUNC__ "lis_matrix_get_diagonal_coo"
lis_matrix_get_diagonal_coo(LIS_MATRIX A,LIS_SCALAR d[])322 LIS_INT lis_matrix_get_diagonal_coo(LIS_MATRIX A, LIS_SCALAR d[])
323 {
324 LIS_INT i;
325 LIS_INT n,nnz;
326
327 LIS_DEBUG_FUNC_IN;
328
329 n = A->n;
330 nnz = A->nnz;
331 if( A->is_splited )
332 {
333 #ifdef _OPENMP
334 #pragma omp parallel for private(i)
335 #endif
336 for(i=0; i<n; i++)
337 {
338 d[i] = A->D->value[i];
339 }
340 }
341 else
342 {
343 for(i=0; i<nnz;i++)
344 {
345 if( A->row[i]==A->col[i] )
346 {
347 d[A->row[i]] = A->value[i];
348 }
349 }
350 }
351 LIS_DEBUG_FUNC_OUT;
352 return LIS_SUCCESS;
353 }
354
355 #undef __FUNC__
356 #define __FUNC__ "lis_matrix_shift_diagonal_coo"
lis_matrix_shift_diagonal_coo(LIS_MATRIX A,LIS_SCALAR sigma)357 LIS_INT lis_matrix_shift_diagonal_coo(LIS_MATRIX A, LIS_SCALAR sigma)
358 {
359 LIS_INT i;
360 LIS_INT n,nnz;
361
362 LIS_DEBUG_FUNC_IN;
363
364 n = A->n;
365 nnz = A->nnz;
366 if( A->is_splited )
367 {
368 #ifdef _OPENMP
369 #pragma omp parallel for private(i)
370 #endif
371 for(i=0; i<n; i++)
372 {
373 A->D->value[i] -= sigma;
374 }
375 }
376 else
377 {
378 for(i=0; i<nnz;i++)
379 {
380 if( A->row[i]==A->col[i] )
381 {
382 A->value[i] -= sigma;
383 }
384 }
385 }
386 LIS_DEBUG_FUNC_OUT;
387 return LIS_SUCCESS;
388 }
389
390 #undef __FUNC__
391 #define __FUNC__ "lis_matrix_scale_coo"
lis_matrix_scale_coo(LIS_MATRIX A,LIS_SCALAR d[])392 LIS_INT lis_matrix_scale_coo(LIS_MATRIX A, LIS_SCALAR d[])
393 {
394 LIS_INT i,j;
395 LIS_INT n,nnz;
396
397 LIS_DEBUG_FUNC_IN;
398
399 n = A->n;
400 if( A->is_splited )
401 {
402 for(j=0; j<A->L->nnz; j++)
403 {
404 i = A->L->row[j];
405 A->L->value[j] *= d[i];
406 }
407 for(i=0;i<n;i++) A->D->value[i] = 1.0;
408 for(j=0; j<A->U->nnz; j++)
409 {
410 i = A->U->row[j];
411 A->U->value[j] *= d[i];
412 }
413 }
414 else
415 {
416 nnz = A->nnz;
417 for(j=0; j<nnz; j++)
418 {
419 i = A->row[j];
420 A->value[j] *= d[i];
421 }
422 }
423 LIS_DEBUG_FUNC_OUT;
424 return LIS_SUCCESS;
425 }
426
427 #undef __FUNC__
428 #define __FUNC__ "lis_matrix_scale_symm_coo"
lis_matrix_scale_symm_coo(LIS_MATRIX A,LIS_SCALAR d[])429 LIS_INT lis_matrix_scale_symm_coo(LIS_MATRIX A, LIS_SCALAR d[])
430 {
431 LIS_INT i,j,k;
432 LIS_INT n,nnz;
433
434 LIS_DEBUG_FUNC_IN;
435
436 n = A->n;
437 if( A->is_splited )
438 {
439 for(k=0; k<A->L->nnz; k++)
440 {
441 i = A->L->row[k];
442 j = A->L->row[k];
443 A->L->value[k] *= d[i]*d[j];
444 }
445 for(i=0;i<n;i++) A->D->value[i] = 1.0;
446 for(k=0; k<A->U->nnz; k++)
447 {
448 i = A->U->row[k];
449 j = A->U->row[k];
450 A->U->value[k] *= d[i]*d[j];
451 }
452 }
453 else
454 {
455 nnz = A->nnz;
456 for(k=0; k<nnz; k++)
457 {
458 i = A->row[k];
459 j = A->row[k];
460 A->value[k] *= d[i]*d[j];
461 }
462 }
463 LIS_DEBUG_FUNC_OUT;
464 return LIS_SUCCESS;
465 }
466
467 #undef __FUNC__
468 #define __FUNC__ "lis_matrix_normf_coo"
lis_matrix_normf_coo(LIS_MATRIX A,LIS_SCALAR * nrm)469 LIS_INT lis_matrix_normf_coo(LIS_MATRIX A, LIS_SCALAR *nrm)
470 {
471 LIS_INT i,j;
472 LIS_INT n;
473 LIS_SCALAR sum;
474
475 LIS_DEBUG_FUNC_IN;
476
477 n = A->n;
478 sum = (LIS_SCALAR)0;
479 if( A->is_splited )
480 {
481 #ifdef _OPENMP
482 #pragma omp parallel for reduction(+:sum) private(i,j)
483 #endif
484 for(i=0; i<n; i++)
485 {
486 sum += A->D->value[i]*A->D->value[i];
487 for(j=A->L->index[i];j<A->L->index[i+1];j++)
488 {
489 sum += A->L->value[j]*A->L->value[j];
490 }
491 for(j=A->U->index[i];j<A->U->index[i+1];j++)
492 {
493 sum += A->U->value[j]*A->U->value[j];
494 }
495 }
496 }
497 else
498 {
499 #ifdef _OPENMP
500 #pragma omp parallel for reduction(+:sum) private(i,j)
501 #endif
502 for(i=0; i<n; i++)
503 {
504 sum += A->value[i]*A->value[i];
505 for(j=A->index[i];j<A->index[i+1];j++)
506 {
507 sum += A->value[j]*A->value[j];
508 }
509 }
510 }
511 *nrm = sqrt(sum);
512 LIS_DEBUG_FUNC_OUT;
513 return LIS_SUCCESS;
514 }
515
516 #undef __FUNC__
517 #define __FUNC__ "lis_matrix_transpose_coo"
lis_matrix_transpose_coo(LIS_MATRIX Ain,LIS_MATRIX * Aout)518 LIS_INT lis_matrix_transpose_coo(LIS_MATRIX Ain, LIS_MATRIX *Aout)
519 {
520
521 LIS_DEBUG_FUNC_IN;
522
523 /* err = lis_matrix_convert_coo2csc(Ain,Aout);*/
524 (*Aout)->matrix_type = LIS_MATRIX_COO;
525 (*Aout)->status = LIS_MATRIX_COO;
526
527 LIS_DEBUG_FUNC_OUT;
528 return LIS_SUCCESS;
529 }
530
531 #undef __FUNC__
532 #define __FUNC__ "lis_matrix_split_coo"
lis_matrix_split_coo(LIS_MATRIX A)533 LIS_INT lis_matrix_split_coo(LIS_MATRIX A)
534 {
535 LIS_INT i,nnz;
536 LIS_INT nnzl,nnzu;
537 LIS_INT err;
538 LIS_INT *lrow,*lcol,*urow,*ucol;
539 LIS_SCALAR *lvalue,*uvalue;
540 LIS_MATRIX_DIAG D;
541
542 LIS_DEBUG_FUNC_IN;
543
544 nnz = A->nnz;
545 nnzl = 0;
546 nnzu = 0;
547 D = NULL;
548 lrow = NULL;
549 lcol = NULL;
550 lvalue = NULL;
551 urow = NULL;
552 ucol = NULL;
553 uvalue = NULL;
554
555 for(i=0;i<nnz;i++)
556 {
557 if( A->col[i]<A->row[i] )
558 {
559 nnzl++;
560 }
561 else if( A->col[i]>A->row[i] )
562 {
563 nnzu++;
564 }
565 }
566
567 err = lis_matrix_LU_create(A);
568 if( err )
569 {
570 return err;
571 }
572 err = lis_matrix_malloc_coo(nnzl,&lrow,&lcol,&lvalue);
573 if( err )
574 {
575 return err;
576 }
577 err = lis_matrix_malloc_coo(nnzu,&urow,&ucol,&uvalue);
578 if( err )
579 {
580 lis_free2(6,lrow,lcol,lvalue,urow,ucol,uvalue);
581 return err;
582 }
583 err = lis_matrix_diag_duplicateM(A,&D);
584 if( err )
585 {
586 lis_free2(6,lrow,lcol,lvalue,urow,ucol,uvalue);
587 return err;
588 }
589
590 nnzl = 0;
591 nnzu = 0;
592 for(i=0;i<nnz;i++)
593 {
594 if( A->col[i]<A->row[i] )
595 {
596 lrow[nnzl] = A->row[i];
597 lcol[nnzl] = A->col[i];
598 lvalue[nnzl] = A->value[i];
599 nnzl++;
600 }
601 else if( A->col[i]>A->row[i] )
602 {
603 urow[nnzu] = A->row[i];
604 ucol[nnzu] = A->col[i];
605 uvalue[nnzu] = A->value[i];
606 nnzu++;
607 }
608 else
609 {
610 D->value[A->row[i]] = A->value[i];
611 }
612 }
613 A->L->nnz = nnzl;
614 A->L->row = lrow;
615 A->L->col = lcol;
616 A->L->value = lvalue;
617 A->U->nnz = nnzu;
618 A->U->row = urow;
619 A->U->col = ucol;
620 A->U->value = uvalue;
621 A->D = D;
622 A->is_splited = LIS_TRUE;
623
624 LIS_DEBUG_FUNC_OUT;
625 return LIS_SUCCESS;
626 }
627
628
629 #undef __FUNC__
630 #define __FUNC__ "lis_matrix_merge_coo"
lis_matrix_merge_coo(LIS_MATRIX A)631 LIS_INT lis_matrix_merge_coo(LIS_MATRIX A)
632 {
633 LIS_INT i;
634 LIS_INT nnz,nnzl,nnzu;
635 LIS_INT err;
636 LIS_INT *row,*col;
637 LIS_SCALAR *value;
638
639 LIS_DEBUG_FUNC_IN;
640
641
642 nnzl = A->L->nnz;
643 nnzu = A->U->nnz;
644 row = NULL;
645 col = NULL;
646 value = NULL;
647 nnz = A->L->nnz + A->U->nnz + A->D->n;
648
649 err = lis_matrix_malloc_coo(nnz,&row,&col,&value);
650 if( err )
651 {
652 return err;
653 }
654
655 nnz = 0;
656 for(i=0;i<nnzl;i++)
657 {
658 row[nnz] = A->L->row[i];
659 col[nnz] = A->L->col[i];
660 value[nnz] = A->L->value[i];
661 nnz++;
662 }
663 for(i=0;i<A->D->n;i++)
664 {
665 row[nnz] = i;
666 col[nnz] = i;
667 value[nnz] = A->D->value[i];
668 nnz++;
669 }
670 for(i=0;i<nnzu;i++)
671 {
672 row[nnz] = A->U->row[i];
673 col[nnz] = A->U->col[i];
674 value[nnz] = A->U->value[i];
675 nnz++;
676 }
677
678 A->nnz = nnz;
679 A->row = row;
680 A->col = col;
681 A->value = value;
682
683 LIS_DEBUG_FUNC_OUT;
684 return LIS_SUCCESS;
685 }
686
687 #undef __FUNC__
688 #define __FUNC__ "lis_matrix_sort_coo"
lis_matrix_sort_coo(LIS_MATRIX A)689 LIS_INT lis_matrix_sort_coo(LIS_MATRIX A)
690 {
691 LIS_INT i,n;
692
693 LIS_DEBUG_FUNC_IN;
694
695 if( !A->is_sorted )
696 {
697 n = A->n;
698 if( A->is_splited )
699 {
700 #ifdef _OPENMP
701 #pragma omp parallel for private(i)
702 #endif
703 for(i=0;i<n;i++)
704 {
705 lis_sort_id(A->L->ptr[i],A->L->ptr[i+1]-1,A->L->index,A->L->value);
706 lis_sort_id(A->U->ptr[i],A->U->ptr[i+1]-1,A->U->index,A->U->value);
707 }
708 }
709 else
710 {
711 #ifdef _OPENMP
712 #pragma omp parallel for private(i)
713 #endif
714 for(i=0;i<n;i++)
715 {
716 lis_sort_id(A->ptr[i],A->ptr[i+1]-1,A->index,A->value);
717 }
718 }
719 A->is_sorted = LIS_TRUE;
720 }
721
722 LIS_DEBUG_FUNC_OUT;
723 return LIS_SUCCESS;
724 }
725
726
727 #undef __FUNC__
728 #define __FUNC__ "lis_matrix_convert_csr2coo"
lis_matrix_convert_csr2coo(LIS_MATRIX Ain,LIS_MATRIX Aout)729 LIS_INT lis_matrix_convert_csr2coo(LIS_MATRIX Ain, LIS_MATRIX Aout)
730 {
731 LIS_INT i,j,k;
732 LIS_INT err;
733 LIS_INT n,nnz;
734 LIS_INT *row,*col;
735 LIS_SCALAR *value;
736
737 LIS_DEBUG_FUNC_IN;
738
739 n = Ain->n;
740 nnz = Ain->nnz;
741
742 row = NULL;
743 col = NULL;
744 value = NULL;
745
746 err = lis_matrix_malloc_coo(nnz,&row,&col,&value);
747 if( err )
748 {
749 return err;
750 }
751
752 /* convert coo */
753 k = 0;
754 for(i=0;i<n;i++)
755 {
756 for(j=Ain->ptr[i];j<Ain->ptr[i+1];j++)
757 {
758 row[k] = i;
759 col[k] = Ain->index[j];
760 value[k] = Ain->value[j];
761 k++;
762 }
763 }
764
765 err = lis_matrix_set_coo(nnz,row,col,value,Aout);
766 if( err )
767 {
768 lis_free2(3,row,col,value);
769 return err;
770 }
771 err = lis_matrix_assemble(Aout);
772 if( err )
773 {
774 lis_matrix_storage_destroy(Aout);
775 return err;
776 }
777
778 LIS_DEBUG_FUNC_OUT;
779 return LIS_SUCCESS;
780 }
781
782 #undef __FUNC__
783 #define __FUNC__ "lis_matrix_convert_coo2csr"
lis_matrix_convert_coo2csr(LIS_MATRIX Ain,LIS_MATRIX Aout)784 LIS_INT lis_matrix_convert_coo2csr(LIS_MATRIX Ain, LIS_MATRIX Aout)
785 {
786 LIS_INT i,j;
787 LIS_INT err;
788 LIS_INT n,nnz;
789 LIS_INT *ptr,*index;
790 LIS_SCALAR *value;
791
792 LIS_DEBUG_FUNC_IN;
793
794 n = Ain->n;
795 nnz = Ain->nnz;
796
797 ptr = NULL;
798 index = NULL;
799 value = NULL;
800
801 err = lis_matrix_malloc_csr(n,nnz,&ptr,&index,&value);
802 if( err )
803 {
804 return err;
805 }
806
807 /* convert csr */
808 lis_sort_iid(0,nnz-1,Ain->row,Ain->col,Ain->value);
809 #ifdef _OPENMP
810 #pragma omp parallel for private(i)
811 #endif
812 for(i=0;i<n+1;i++)
813 {
814 ptr[i] = 0;
815 }
816 for(i=0;i<nnz;i++)
817 {
818 ptr[Ain->row[i]+1]++;
819 }
820 for(i=0;i<n;i++)
821 {
822 ptr[i+1] += ptr[i];
823 }
824 #ifdef _OPENMP
825 #pragma omp parallel for private(i,j)
826 #endif
827 for(i=0;i<n;i++)
828 {
829 for(j=ptr[i];j<ptr[i+1];j++)
830 {
831 value[j] = Ain->value[j];
832 index[j] = Ain->col[j];
833 }
834 }
835
836 err = lis_matrix_set_csr(nnz,ptr,index,value,Aout);
837 if( err )
838 {
839 lis_free2(3,ptr,index,value);
840 return err;
841 }
842 err = lis_matrix_assemble(Aout);
843 if( err )
844 {
845 lis_matrix_storage_destroy(Aout);
846 return err;
847 }
848
849 LIS_DEBUG_FUNC_OUT;
850 return LIS_SUCCESS;
851 }
852