Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
BlaFormat.c
Go to the documentation of this file.
1
15#include "fasp.h"
16#include "fasp_block.h"
17#include "fasp_functs.h"
18
19/*---------------------------------*/
20/*-- Public Functions --*/
21/*---------------------------------*/
22
36SHORT fasp_format_dcoo_dcsr (const dCOOmat *A,
37 dCSRmat *B)
38{
39 const INT m=A->row, n=A->col, nnz=A->nnz;
40 INT iind, jind, i;
41
42 fasp_dcsr_alloc(m,n,nnz,B);
43 INT *ia = B->IA;
44
45 INT *ind = (INT *) fasp_mem_calloc(m+1,sizeof(INT));
46 memset(ind, 0, sizeof(INT)*(m+1)); // initialize ind
47 for ( i=0; i<nnz; ++i ) ind[A->rowind[i]+1]++; // count nnz in each row
48
49 ia[0] = 0; // first index starting from zero
50 for ( i=1; i<=m; ++i ) {
51 ia[i] = ia[i-1]+ind[i]; // set row_idx
52 ind[i] = ia[i];
53 }
54
55 // loop over nnz and set col_idx and val
56 for ( i=0; i<nnz; ++i ) {
57 iind = A->rowind[i]; jind = ind[iind];
58 B->JA [jind] = A->colind[i];
59 B->val[jind] = A->val[i];
60 ind[iind] = ++jind;
61 }
62
63 fasp_mem_free(ind); ind = NULL;
64
65 return FASP_SUCCESS;
66}
67
83SHORT fasp_format_dcsr_dcoo (const dCSRmat *A,
84 dCOOmat *B)
85{
86 const INT m=A->row, nnz=A->nnz;
87 INT i, j;
88
89 B->rowind = (INT *)fasp_mem_calloc(nnz,sizeof(INT));
90 B->colind = (INT *)fasp_mem_calloc(nnz,sizeof(INT));
91 B->val = (REAL *)fasp_mem_calloc(nnz,sizeof(REAL));
92
93#ifdef _OPENMP
94#pragma omp parallel for if(m>OPENMP_HOLDS) private(i, j)
95#endif
96 for (i=0;i<m;++i) {
97 for (j=A->IA[i];j<A->IA[i+1];++j) B->rowind[j]=i;
98 }
99
100 memcpy(B->colind, A->JA, nnz*sizeof(INT));
101 memcpy(B->val, A->val, nnz*sizeof(REAL));
102
103 return FASP_SUCCESS;
104}
105
119SHORT fasp_format_dstr_dcsr (const dSTRmat *A,
120 dCSRmat *B)
121{
122 // some members of A
123 const INT nc = A->nc;
124 const INT ngrid = A->ngrid;
125 const INT nband = A->nband;
126 const INT *offsets = A->offsets;
127
128 REAL *diag = A->diag;
129 REAL **offdiag = A->offdiag;
130
131 // some members of B
132 const INT glo_row = nc*ngrid;
133 INT glo_nnz;
134 INT *ia = NULL;
135 INT *ja = NULL;
136 REAL *a = NULL;
137
138 dCSRmat B_tmp;
139
140 // local variables
141 INT width;
142 INT nc2 = nc*nc;
143 INT BAND,ROW,COL;
144 INT ncb,nci;
145 INT row_start,col_start;
146 INT block; // how many blocks in the current ROW
147 INT i,j;
148 INT pos;
149 INT start;
150 INT val_L_start,val_R_start;
151 INT row;
152 INT tmp_col;
153 REAL tmp_val;
154
155 // allocate for 'ia' array
156 ia = (INT *)fasp_mem_calloc(glo_row+1,sizeof(INT));
157
158 // Generate the 'ia' array
159 ia[0] = 0;
160 for (ROW = 0; ROW < ngrid; ++ROW) {
161 block = 1; // diagonal block
162 for (BAND = 0; BAND < nband; ++BAND) {
163 width = offsets[BAND];
164 COL = ROW + width;
165 if (width < 0) {
166 if (COL >= 0) ++block;
167 }
168 else {
169 if (COL < ngrid) ++block;
170 }
171 } // end for BAND
172
173 ncb = nc*block;
174 row_start = ROW*nc;
175
176 for (i = 0; i < nc; i ++) {
177 row = row_start + i;
178 ia[row+1] = ia[row] + ncb;
179 }
180 } // end for ROW
181
182 // allocate for 'ja' and 'a' arrays
183 glo_nnz = ia[glo_row];
184 ja = (INT *)fasp_mem_calloc(glo_nnz,sizeof(INT));
185 a = (REAL *)fasp_mem_calloc(glo_nnz,sizeof(REAL));
186
187 // Generate the 'ja' and 'a' arrays at the same time
188 for (ROW = 0; ROW < ngrid; ++ROW) {
189 row_start = ROW*nc;
190 val_L_start = ROW*nc2;
191
192 // deal with the diagonal band
193 for (i = 0; i < nc; i ++) {
194 nci = nc*i;
195 row = row_start + i;
196 start = ia[row];
197 for (j = 0; j < nc; j ++) {
198 pos = start + j;
199 ja[pos] = row_start + j;
200 a[pos] = diag[val_L_start+nci+j];
201 }
202 }
203 block = 1;
204
205 // deal with the off-diagonal bands
206 for (BAND = 0; BAND < nband; ++BAND) {
207 width = offsets[BAND];
208 COL = ROW + width;
209 ncb = nc*block;
210 col_start = COL*nc;
211
212 if (width < 0) {
213 if (COL >= 0) {
214 val_R_start = COL*nc2;
215 for (i = 0; i < nc; i ++) {
216 nci = nc*i;
217 row = row_start + i;
218 start = ia[row];
219 for (j = 0 ; j < nc; j ++) {
220 pos = start + ncb + j;
221 ja[pos] = col_start + j;
222 a[pos] = offdiag[BAND][val_R_start+nci+j];
223 }
224 }
225 ++block;
226 }
227 }
228 else {
229 if (COL < ngrid) {
230 for (i = 0; i < nc; i ++) {
231 nci = nc*i;
232 row = row_start + i;
233 start = ia[row];
234 for (j = 0; j < nc; j ++) {
235 pos = start + ncb + j;
236 ja[pos] = col_start + j;
237 a[pos] = offdiag[BAND][val_L_start+nci+j];
238 }
239 }
240 ++block;
241 }
242 }
243 }
244 }
245
246 // Reordering in such manner that every diagonal element
247 // is firstly stored in the corresponding row
248 if (nc > 1) {
249 for (ROW = 0; ROW < ngrid; ++ROW) {
250 row_start = ROW*nc;
251 for (j = 1; j < nc; j ++) {
252 row = row_start + j;
253 start = ia[row];
254 pos = start + j;
255
256 // swap in 'ja'
257 tmp_col = ja[start];
258 ja[start] = ja[pos];
259 ja[pos] = tmp_col;
260
261 // swap in 'a'
262 tmp_val = a[start];
263 a[start] = a[pos];
264 a[pos] = tmp_val;
265 }
266 }
267 }
268
269 /* fill all the members of B_tmp */
270 B_tmp.row = glo_row;
271 B_tmp.col = glo_row;
272 B_tmp.nnz = glo_nnz;
273 B_tmp.IA = ia;
274 B_tmp.JA = ja;
275 B_tmp.val = a;
276
277 *B = B_tmp;
278
279 return FASP_SUCCESS;
280}
281
294dCSRmat fasp_format_dblc_dcsr (const dBLCmat *Ab)
295{
296 const INT mb=Ab->brow, nb=Ab->bcol, nbl=mb*nb;
297 dCSRmat **blockptr=Ab->blocks, *blockptrij, A;
298
299 INT i,j,ij,ir,i1,length,ilength,start,irmrow,irmrow1;
300 INT *row, *col;
301 INT m=0,n=0,nnz=0;
302
303 row = (INT *)fasp_mem_calloc(mb+1,sizeof(INT));
304 col = (INT *)fasp_mem_calloc(nb+1,sizeof(INT));
305
306 // count the size of A
307 row[0]=0; col[0]=0;
308 for (i=0;i<mb;++i) { m+=blockptr[i*nb]->row; row[i+1]=m; }
309 for (i=0;i<nb;++i) { n+=blockptr[i]->col; col[i+1]=n; }
310
311#ifdef _OPENMP
312#pragma omp parallel for reduction(+:nnz) if (nbl>OPENMP_HOLDS) private(i)
313#endif
314 for (i=0;i<nbl;++i) { nnz+=blockptr[i]->nnz; }
315
316 // memory space allocation
317 A = fasp_dcsr_create(m,n,nnz);
318
319 // set dCSRmat for A
320 A.IA[0]=0;
321 for (i=0;i<mb;++i) {
322
323 for (ir=row[i];ir<row[i+1];ir++) {
324
325 for (length=j=0;j<nb;++j) {
326 ij=i*nb+j; blockptrij=blockptr[ij];
327 if (blockptrij->nnz>0) {
328 start=A.IA[ir]+length;
329 irmrow=ir-row[i]; irmrow1=irmrow+1;
330 ilength=blockptrij->IA[irmrow1]-blockptrij->IA[irmrow];
331 if (ilength>0) {
332 memcpy(&(A.val[start]),&(blockptrij->val[blockptrij->IA[irmrow]]),ilength*sizeof(REAL));
333 memcpy(&(A.JA[start]), &(blockptrij->JA[blockptrij->IA[irmrow]]), ilength*sizeof(INT));
334 for (i1=0;i1<ilength;i1++) A.JA[start+i1]+=col[j];
335 length+=ilength;
336 }
337 }
338 } // end for j
339
340 A.IA[ir+1]=A.IA[ir]+length;
341 } // end for ir
342
343 } // end for i
344
345 fasp_mem_free(row); row = NULL;
346 fasp_mem_free(col); col = NULL;
347
348 return(A);
349}
350
363dCSRLmat * fasp_format_dcsrl_dcsr (const dCSRmat *A)
364{
365 REAL *DATA = A -> val;
366 INT *IA = A -> IA;
367 INT *JA = A -> JA;
368 INT num_rows = A -> row;
369 INT num_cols = A -> col;
370 INT num_nonzeros = A -> nnz;
371
372 dCSRLmat *B = NULL;
373 INT dif;
374 INT *nzdifnum = NULL;
375 INT *rowstart = NULL;
376 INT *rowindex = (INT *)fasp_mem_calloc(num_rows, sizeof(INT));
377 INT *ja = (INT *)fasp_mem_calloc(num_nonzeros, sizeof(INT));
378 REAL *data = (REAL *)fasp_mem_calloc(num_nonzeros, sizeof(REAL));
379
380 /* auxiliary arrays */
381 INT *nzrow = (INT *)fasp_mem_calloc(num_rows, sizeof(INT));
382 INT *counter = NULL;
383 INT *invnzdif = NULL;
384
385 INT i,j,k,cnt,maxnzrow;
386
387 //-----------------------------------------
388 // Generate 'nzrow' and 'maxnzrow'
389 //-----------------------------------------
390
391 maxnzrow = 0;
392 for (i = 0; i < num_rows; i ++) {
393 nzrow[i] = IA[i+1] - IA[i];
394 if (nzrow[i] > maxnzrow) {
395 maxnzrow = nzrow[i];
396 }
397 }
398 /* generate 'counter' */
399 counter = (INT *)fasp_mem_calloc(maxnzrow + 1, sizeof(INT));
400
401 for (i = 0; i < num_rows; i ++) {
402 counter[nzrow[i]] ++;
403 }
404
405 //--------------------------------------------
406 // Determine 'dif'
407 //--------------------------------------------
408
409 for (dif = 0, i = 0; i < maxnzrow + 1; i ++) {
410 if (counter[i] > 0) dif ++;
411 }
412
413 //--------------------------------------------
414 // Generate the 'nzdifnum' and 'rowstart'
415 //--------------------------------------------
416
417 nzdifnum = (INT *)fasp_mem_calloc(dif, sizeof(INT));
418 invnzdif = (INT *)fasp_mem_calloc(maxnzrow + 1, sizeof(INT));
419 rowstart = (INT *)fasp_mem_calloc(dif + 1, sizeof(INT));
420 rowstart[0] = 0;
421 for (cnt = 0, i = 0; i < maxnzrow + 1; i ++) {
422 if (counter[i] > 0) {
423 nzdifnum[cnt] = i;
424 invnzdif[i] = cnt;
425 rowstart[cnt+1] = rowstart[cnt] + counter[i];
426 cnt ++;
427 }
428 }
429
430 //--------------------------------------------
431 // Generate the 'rowindex'
432 //--------------------------------------------
433
434 for (i = 0; i < num_rows; i ++) {
435 j = invnzdif[nzrow[i]];
436 rowindex[rowstart[j]] = i;
437 rowstart[j] ++;
438 }
439 /* recover 'rowstart' */
440 for (i = dif; i > 0; i --) {
441 rowstart[i] = rowstart[i-1];
442 }
443 rowstart[0] = 0;
444
445 //--------------------------------------------
446 // Generate the 'data' and 'ja'
447 //--------------------------------------------
448
449 for (cnt = 0, i = 0; i < num_rows; i ++) {
450 k = rowindex[i];
451 for (j = IA[k]; j < IA[k+1]; j ++) {
452 data[cnt] = DATA[j];
453 ja[cnt] = JA[j];
454 cnt ++;
455 }
456 }
457
458 //------------------------------------------------------------
459 // Create and fill a dCSRLmat B
460 //------------------------------------------------------------
461
462 B = fasp_dcsrl_create(num_rows, num_cols, num_nonzeros);
463 B -> dif = dif;
464 B -> nz_diff = nzdifnum;
465 B -> index = rowindex;
466 B -> start = rowstart;
467 B -> ja = ja;
468 B -> val = data;
469
470 //----------------------------
471 // Free the auxiliary arrays
472 //----------------------------
473
474 free(nzrow);
475 free(counter);
476 free(invnzdif);
477
478 return B;
479}
480
497dCSRmat fasp_format_dbsr_dcsr (const dBSRmat *B)
498{
499 dCSRmat A;
500
501 /* members of B */
502 INT ROW = B->ROW;
503 INT COL = B->COL;
504 INT NNZ = B->NNZ;
505 INT nb = B->nb;
506 INT *IA = B->IA;
507 INT *JA = B->JA;
508 REAL *val = B->val;
509
510 INT storage_manner = B->storage_manner;
511
512 INT jump = nb*nb;
513 INT rowA = ROW*nb;
514 INT colA = COL*nb;
515 INT nzA = NNZ*jump;
516
517 INT *ia = NULL;
518 INT *ja = NULL;
519 REAL *a = NULL;
520
521 INT i,j,k;
522 INT mr,mc;
523 INT rowstart0,rowstart,colstart0,colstart;
524 INT colblock,nzperrow;
525
526 REAL *vp = NULL;
527 REAL *ap = NULL;
528 INT *jap = NULL;
529
530 SHORT use_openmp = FALSE;
531
532#ifdef _OPENMP
533 INT stride_i,mybegin,myend,myid,nthreads;
534 if ( ROW > OPENMP_HOLDS ) {
535 use_openmp = TRUE;
536 nthreads = fasp_get_num_threads();
537 }
538#endif
539
540 //--------------------------------------------------------
541 // Create a CSR Matrix
542 //--------------------------------------------------------
543 A = fasp_dcsr_create(rowA, colA, nzA);
544 ia = A.IA;
545 ja = A.JA;
546 a = A.val;
547
548 //--------------------------------------------------------------------------
549 // Compute the number of nonzeros per row, and after this loop,
550 // ia[i],i=1:rowA, will be the number of nonzeros of the (i-1)-th row.
551 //--------------------------------------------------------------------------
552
553 if (use_openmp) {
554#ifdef _OPENMP
555 stride_i = ROW/nthreads;
556#pragma omp parallel private(myid, mybegin, myend, i, rowstart, colblock, nzperrow, j)
557 {
558 myid = omp_get_thread_num();
559 mybegin = myid*stride_i;
560 if(myid < nthreads-1) myend = mybegin+stride_i;
561 else myend = ROW;
562 for (i=mybegin; i<myend; ++i)
563 {
564 rowstart = i*nb + 1;
565 colblock = IA[i+1] - IA[i];
566 nzperrow = colblock*nb;
567 for (j = 0; j < nb; ++j)
568 {
569 ia[rowstart+j] = nzperrow;
570 }
571 }
572 }
573#endif
574 }
575 else {
576 for (i = 0; i < ROW; ++i)
577 {
578 rowstart = i*nb + 1;
579 colblock = IA[i+1] - IA[i];
580 nzperrow = colblock*nb;
581 for (j = 0; j < nb; ++j)
582 {
583 ia[rowstart+j] = nzperrow;
584 }
585 }
586 }
587
588 //-----------------------------------------------------
589 // Generate the real 'ia' for CSR of A
590 //-----------------------------------------------------
591
592 ia[0] = 0;
593 for (i = 1; i <= rowA; ++i)
594 {
595 ia[i] += ia[i-1];
596 }
597
598 //-----------------------------------------------------
599 // Generate 'ja' and 'a' for CSR of A
600 //-----------------------------------------------------
601
602 switch (storage_manner)
603 {
604 case 0: // each non-zero block elements are stored in row-major order
605 {
606 if (use_openmp) {
607#ifdef _OPENMP
608#pragma omp parallel private(myid, mybegin, myend, i, k, j, rowstart, colstart, vp, mr, ap, jap, mc)
609 {
610 myid = omp_get_thread_num();
611 mybegin = myid*stride_i;
612 if(myid < nthreads-1) myend = mybegin+stride_i;
613 else myend = ROW;
614 for (i=mybegin; i<myend; ++i)
615 {
616 for (k = IA[i]; k < IA[i+1]; ++k)
617 {
618 j = JA[k];
619 rowstart = i*nb;
620 colstart = j*nb;
621 vp = &val[k*jump];
622 for (mr = 0; mr < nb; mr ++)
623 {
624 ap = &a[ia[rowstart]];
625 jap = &ja[ia[rowstart]];
626 for (mc = 0; mc < nb; mc ++)
627 {
628 *ap = *vp;
629 *jap = colstart + mc;
630 vp ++; ap ++; jap ++;
631 }
632 ia[rowstart] += nb;
633 rowstart ++;
634 }
635 }
636 }
637 }
638#endif
639 }
640 else {
641 for (i = 0; i < ROW; ++i)
642 {
643 for (k = IA[i]; k < IA[i+1]; ++k)
644 {
645 j = JA[k];
646 rowstart = i*nb;
647 colstart = j*nb;
648 vp = &val[k*jump];
649 for (mr = 0; mr < nb; mr ++)
650 {
651 ap = &a[ia[rowstart]];
652 jap = &ja[ia[rowstart]];
653 for (mc = 0; mc < nb; mc ++)
654 {
655 *ap = *vp;
656 *jap = colstart + mc;
657 vp ++; ap ++; jap ++;
658 }
659 ia[rowstart] += nb;
660 rowstart ++;
661 }
662 }
663 }
664 }
665 }
666 break;
667
668 case 1: // each non-zero block elements are stored in column-major order
669 {
670 for (i = 0; i < ROW; ++i)
671 {
672 for (k = IA[i]; k < IA[i+1]; ++k)
673 {
674 j = JA[k];
675 rowstart0 = i*nb;
676 colstart0 = j*nb;
677 vp = &val[k*jump];
678 for (mc = 0; mc < nb; mc ++)
679 {
680 rowstart = rowstart0;
681 colstart = colstart0 + mc;
682 for (mr = 0; mr < nb; mr ++)
683 {
684 a[ia[rowstart]] = *vp;
685 ja[ia[rowstart]] = colstart;
686 vp ++; ia[rowstart]++; rowstart++;
687 }
688 }
689 }
690 }
691 }
692 break;
693 }
694
695 //-----------------------------------------------------
696 // Map back the real 'ia' for CSR of A
697 //-----------------------------------------------------
698
699 for (i = rowA; i > 0; i --) {
700 ia[i] = ia[i-1];
701 }
702 ia[0] = 0;
703
704 return (A);
705}
706
723dBSRmat fasp_format_dcsr_dbsr (const dCSRmat *A,
724 const INT nb)
725{
726 INT i, j, k, ii, jj, kk, l, mod, nnz;
727 INT row = A->row/nb;
728 INT col = A->col/nb;
729 INT nb2 = nb*nb;
730 INT *IA = A->IA;
731 INT *JA = A->JA;
732 REAL *val = A->val;
733
734 dBSRmat B; // Safe-guard check
735 INT *col_flag, *ia, *ja;
736 REAL *bval;
737
738 if ((A->row)%nb!=0) {
739 printf("### ERROR: A.row=%d is not a multiplication of nb=%d!\n",
740 A->row, nb);
741 fasp_chkerr(ERROR_MAT_SIZE, __FUNCTION__);
742 }
743
744 if ((A->col)%nb!=0) {
745 printf("### ERROR: A.col=%d is not a multiplication of nb=%d!\n",
746 A->col, nb);
747 fasp_chkerr(ERROR_MAT_SIZE, __FUNCTION__);
748 }
749
750 B.ROW = row;
751 B.COL = col;
752 B.nb = nb;
753 B.storage_manner = 0;
754
755 // allocate memory for B
756 col_flag = (INT *)fasp_mem_calloc(col, sizeof(INT));
757 ia = (INT *) fasp_mem_calloc(row+1, sizeof(INT));
758
759 fasp_iarray_set(col, col_flag, -1);
760
761 // Get ia for BSR format
762 nnz = 0;
763 for (i=0; i<row; ++i) {
764 ii = nb*i;
765 for (j=0; j<nb; ++j) {
766 jj = ii+j;
767 for (k=IA[jj]; k<IA[jj+1]; ++k) {
768 kk = JA[k]/nb;
769 if (col_flag[kk]!=0) {
770 col_flag[kk] = 0;
771 //ja[nnz] = kk;
772 nnz ++;
773 }
774 }
775 }
776 ia[i+1] = nnz;
777 fasp_iarray_set(col, col_flag, -1);
778 }
779
780 // set NNZ
781 B.NNZ = nnz;
782
783 // allocate ja and bval
784 ja = (INT*)fasp_mem_calloc(nnz, sizeof(INT));
785 bval = (REAL*)fasp_mem_calloc(nnz*nb2, sizeof(REAL));
786
787 // Get ja for BSR format
788 nnz = 0;
789 for (i=0; i<row; ++i) {
790 ii = nb*i;
791 for(j=0; j<nb; ++j) {
792 jj = ii+j;
793 for(k=IA[jj]; k<IA[jj+1]; ++k) {
794 kk = JA[k]/nb;
795 if (col_flag[kk]!=0) {
796 col_flag[kk] = 0;
797 ja[nnz] = kk;
798 nnz ++;
799 }
800 }
801 }
802 ia[i+1] = nnz;
803 fasp_iarray_set(col, col_flag, -1);
804 }
805
806 // Get non-zeros of BSR
807 for (i=0; i<row; ++i) {
808 ii = nb*i;
809 for (j=0; j<nb; ++j) {
810 jj = ii+j;
811 for (k=IA[jj]; k<IA[jj+1]; ++k) {
812 for (l=ia[i]; l<ia[i+1]; ++l) {
813 if (JA[k]/nb ==ja[l]) {
814 mod = JA[k]%nb;
815 bval[l*nb2+j*nb+mod] = val[k];
816 break;
817 }
818 }
819 }
820 }
821 }
822
823 B.IA = ia;
824 B.JA = ja;
825 B.val = bval;
826
827 fasp_mem_free(col_flag); col_flag = NULL;
828
829 return B;
830}
831
844dBSRmat fasp_format_dstr_dbsr (const dSTRmat *B)
845{
846 // members of 'B'
847 INT nc = B->nc;
848 INT ngrid = B->ngrid;
849 REAL *diag = B->diag;
850 INT nband = B->nband;
851 INT *offsets = B->offsets;
852 REAL **offdiag = B->offdiag;
853
854 // members of 'A'
855 dBSRmat A;
856 INT NNZ;
857 INT *IA = NULL;
858 INT *JA = NULL;
859 REAL *val = NULL;
860
861 // local variables
862 INT i,j,k,m;
863 INT nc2 = nc*nc;
864 INT ngridplus1 = ngrid + 1;
865
866 // compute NNZ
867 NNZ = ngrid;
868 for (i = 0; i < nband; ++i) {
869 NNZ += (ngrid - abs(offsets[i]));
870 }
871
872 // Create and Initialize a dBSRmat 'A'
873 A = fasp_dbsr_create(ngrid, ngrid, NNZ, nc, 0);
874 IA = A.IA;
875 JA = A.JA;
876 val = A.val;
877
878 // Generate 'IA'
879 for (i = 1; i < ngridplus1; ++i) IA[i] = 1; // take the diagonal blocks into account
880 for (i = 0; i < nband; ++i) {
881 k = offsets[i];
882 if (k < 0) {
883 for (j = -k+1; j < ngridplus1; ++j) {
884 IA[j] ++;
885 }
886 }
887 else {
888 m = ngridplus1 - k;
889 for (j = 1; j < m; ++j)
890 {
891 IA[j] ++;
892 }
893 }
894 }
895 IA[0] = 0;
896 for (i = 1; i < ngridplus1; ++i) {
897 IA[i] += IA[i-1];
898 }
899
900 // Generate 'JA' and 'val' at the same time
901 for (i = 0 ; i < ngrid; ++i) {
902 memcpy(val + IA[i]*nc2, diag + i*nc2, nc2*sizeof(REAL));
903 JA[IA[i]] = i;
904 IA[i] ++;
905 }
906
907 for (i = 0; i < nband; ++i) {
908 k = offsets[i];
909 if (k < 0) {
910 for (j = -k; j < ngrid; ++j) {
911 m = j + k;
912 memcpy(val+IA[j]*nc2, offdiag[i]+m*nc2, nc2*sizeof(REAL));
913 JA[IA[j]] = m;
914 IA[j] ++;
915 }
916 }
917 else {
918 m = ngrid - k;
919 for (j = 0; j < m; ++j) {
920 memcpy(val + IA[j]*nc2, offdiag[i] + j*nc2, nc2*sizeof(REAL));
921 JA[IA[j]] = k + j;
922 IA[j] ++;
923 }
924 }
925 }
926
927 // Map back the real 'IA' for BSR of A
928 for (i = ngrid; i > 0; i --) {
929 IA[i] = IA[i-1];
930 }
931 IA[0] = 0;
932
933 return (A);
934}
935
948dCOOmat * fasp_format_dbsr_dcoo (const dBSRmat *B)
949{
950 /* members of B */
951 INT ROW = B->ROW;
952 INT COL = B->COL;
953 INT NNZ = B->NNZ;
954 INT nb = B->nb;
955 INT *IA = B->IA;
956 INT *JA = B->JA;
957 REAL *val = B->val;
958
959 dCOOmat *A = NULL;
960 INT nb2 = nb*nb;
961 INT num_nonzeros = NNZ*nb2;
962 INT *rowA = NULL;
963 INT *colA = NULL;
964 REAL *valA = NULL;
965
966 INT i,j,k,inb;
967 INT row_start, col_start;
968 INT cnt,mr,mc;
969 REAL *pt = NULL;
970
971 // Create and Initialize a dCOOmat 'A'
972 A = (dCOOmat *)fasp_mem_calloc(1, sizeof(dCOOmat));
973 A->row = ROW*nb;
974 A->col = COL*nb;
975 A->nnz = num_nonzeros;
976 rowA = (INT *)fasp_mem_calloc(num_nonzeros, sizeof(INT));
977 colA = (INT *)fasp_mem_calloc(num_nonzeros, sizeof(INT));
978 valA = (REAL *)fasp_mem_calloc(num_nonzeros, sizeof(REAL));
979 A->rowind = rowA;
980 A->colind = colA;
981 A->val = valA;
982
983 cnt = 0;
984 for (i = 0; i < ROW; ++i) {
985 inb = i*nb;
986 for (k = IA[i]; k < IA[i+1]; ++k) {
987 j = JA[k];
988 pt = &val[k*nb2];
989 row_start = inb;
990 col_start = j*nb;
991 for (mr = 0; mr < nb; mr ++) {
992 for (mc = 0; mc < nb; mc ++) {
993 rowA[cnt] = row_start;
994 colA[cnt] = col_start + mc;
995 valA[cnt] = (*pt);
996 pt ++;
997 cnt ++;
998 }
999 row_start ++;
1000 }
1001 }
1002 }
1003
1004 return (A);
1005}
1006
1007/*---------------------------------*/
1008/*-- End of File --*/
1009/*---------------------------------*/
fasp_iarray_set
void fasp_iarray_set(const INT n, INT *x, const INT val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:98
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_chkerr
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
fasp_format_dcsr_dbsr
dBSRmat fasp_format_dcsr_dbsr(const dCSRmat *A, const INT nb)
Transfer a dCSRmat type matrix into a dBSRmat.
Definition: BlaFormat.c:723
fasp_format_dcsr_dcoo
SHORT fasp_format_dcsr_dcoo(const dCSRmat *A, dCOOmat *B)
Transform a REAL matrix from its CSR format to its IJ format.
Definition: BlaFormat.c:83
fasp_format_dstr_dcsr
SHORT fasp_format_dstr_dcsr(const dSTRmat *A, dCSRmat *B)
Transfer a 'dSTRmat' type matrix into a 'dCSRmat' type matrix.
Definition: BlaFormat.c:119
fasp_format_dcsrl_dcsr
dCSRLmat * fasp_format_dcsrl_dcsr(const dCSRmat *A)
Convert a dCSRmat into a dCSRLmat.
Definition: BlaFormat.c:363
fasp_format_dblc_dcsr
dCSRmat fasp_format_dblc_dcsr(const dBLCmat *Ab)
Form the whole dCSRmat A using blocks given in Ab.
Definition: BlaFormat.c:294
fasp_format_dcoo_dcsr
SHORT fasp_format_dcoo_dcsr(const dCOOmat *A, dCSRmat *B)
Transform a REAL matrix from its IJ format to its CSR format.
Definition: BlaFormat.c:36
fasp_format_dbsr_dcoo
dCOOmat * fasp_format_dbsr_dcoo(const dBSRmat *B)
Transfer a 'dBSRmat' type matrix to a 'dCOOmat' type matrix.
Definition: BlaFormat.c:948
fasp_format_dbsr_dcsr
dCSRmat fasp_format_dbsr_dcsr(const dBSRmat *B)
Transfer a 'dBSRmat' type matrix into a dCSRmat.
Definition: BlaFormat.c:497
fasp_format_dstr_dbsr
dBSRmat fasp_format_dstr_dbsr(const dSTRmat *B)
Transfer a 'dSTRmat' type matrix to a 'dBSRmat' type matrix.
Definition: BlaFormat.c:844
ilength
int ilength
Definition: BlaIO.c:23
fasp_dbsr_create
dBSRmat fasp_dbsr_create(const INT ROW, const INT COL, const INT NNZ, const INT nb, const INT storage_manner)
Create BSR sparse matrix data memory space.
Definition: BlaSparseBSR.c:48
fasp_dcsrl_create
dCSRLmat * fasp_dcsrl_create(const INT num_rows, const INT num_cols, const INT num_nonzeros)
Create a dCSRLmat object.
Definition: BlaSparseCSRL.c:39
fasp_dcsr_create
dCSRmat fasp_dcsr_create(const INT m, const INT n, const INT nnz)
Create CSR sparse matrix data memory space.
Definition: BlaSparseCSR.c:47
fasp_dcsr_alloc
void fasp_dcsr_alloc(const INT m, const INT n, const INT nnz, dCSRmat *A)
Allocate CSR sparse matrix memory space.
Definition: BlaSparseCSR.c:138
fasp.h
Main header file for the FASP project.
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
INT
#define INT
Definition: fasp.h:72
fasp_block.h
Header file for FASP block matrices.
FASP_SUCCESS
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:19
OPENMP_HOLDS
#define OPENMP_HOLDS
Definition: fasp_const.h:269
TRUE
#define TRUE
Definition of logic type.
Definition: fasp_const.h:61
FALSE
#define FALSE
Definition: fasp_const.h:62
ERROR_MAT_SIZE
#define ERROR_MAT_SIZE
Definition: fasp_const.h:26
dBLCmat
Block REAL CSR matrix format.
Definition: fasp_block.h:74
dBLCmat::brow
INT brow
row number of blocks in A, m
Definition: fasp_block.h:77
dBLCmat::blocks
dCSRmat ** blocks
blocks of dCSRmat, point to blocks[brow][bcol]
Definition: fasp_block.h:83
dBLCmat::bcol
INT bcol
column number of blocks A, n
Definition: fasp_block.h:80
dBSRmat
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:34
dBSRmat::COL
INT COL
number of cols of sub-blocks in matrix A, N
Definition: fasp_block.h:40
dBSRmat::NNZ
INT NNZ
number of nonzero sub-blocks in matrix A, NNZ
Definition: fasp_block.h:43
dBSRmat::val
REAL * val
Definition: fasp_block.h:57
dBSRmat::nb
INT nb
dimension of each sub-block
Definition: fasp_block.h:46
dBSRmat::IA
INT * IA
integer array of row pointers, the size is ROW+1
Definition: fasp_block.h:60
dBSRmat::ROW
INT ROW
number of rows of sub-blocks in matrix A, M
Definition: fasp_block.h:37
dBSRmat::JA
INT * JA
Definition: fasp_block.h:64
dBSRmat::storage_manner
INT storage_manner
storage manner for each sub-block
Definition: fasp_block.h:49
dCOOmat
Sparse matrix of REAL type in COO (IJ) format.
Definition: fasp.h:221
dCOOmat::colind
INT * colind
integer array of column indices, the size is nnz
Definition: fasp.h:236
dCOOmat::col
INT col
column of matrix A, n
Definition: fasp.h:227
dCOOmat::rowind
INT * rowind
integer array of row indices, the size is nnz
Definition: fasp.h:233
dCOOmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:239
dCOOmat::row
INT row
row number of matrix A, m
Definition: fasp.h:224
dCOOmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:230
dCSRLmat
Sparse matrix of REAL type in CSRL format.
Definition: fasp.h:277
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::col
INT col
column of matrix A, n
Definition: fasp.h:157
dCSRmat::val
REAL * val
nonzero entries of A
Definition: fasp.h:169
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dCSRmat::IA
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:163
dCSRmat::nnz
INT nnz
number of nonzero entries
Definition: fasp.h:160
dCSRmat::JA
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:166
dSTRmat
Structure matrix of REAL type.
Definition: fasp.h:316
dSTRmat::offsets
INT * offsets
offsets of the off-diagonals (length is nband)
Definition: fasp.h:343
dSTRmat::nband
INT nband
number of off-diag bands
Definition: fasp.h:340
dSTRmat::diag
REAL * diag
diagonal entries (length is ngrid*(nc^2))
Definition: fasp.h:337
dSTRmat::ngrid
INT ngrid
number of grids
Definition: fasp.h:334
dSTRmat::nc
INT nc
size of each block (number of components)
Definition: fasp.h:331
dSTRmat::offdiag
REAL ** offdiag
off-diagonal entries (dimension is nband * [(ngrid-|offsets|) * nc^2])
Definition: fasp.h:346