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