Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
sparse_csr.c
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <time.h>
8 
9 #ifdef _OPENMP
10 #include <omp.h>
11 #endif
12 
13 #include "fasp.h"
14 #include "fasp_functs.h"
15 
16 /*---------------------------------*/
17 /*-- Public Functions --*/
18 /*---------------------------------*/
19 
35  const INT n,
36  const INT nnz)
37 {
38  dCSRmat A;
39 
40  if ( m > 0 ) {
41  A.IA = (INT *)fasp_mem_calloc(m+1, sizeof(INT));
42  }
43  else {
44  A.IA = NULL;
45  }
46 
47  if ( n > 0 ) {
48  A.JA = (INT *)fasp_mem_calloc(nnz, sizeof(INT));
49  }
50  else {
51  A.JA = NULL;
52  }
53 
54  if ( nnz > 0 ) {
55  A.val = (REAL *)fasp_mem_calloc(nnz, sizeof(REAL));
56  }
57  else {
58  A.val = NULL;
59  }
60 
61  A.row = m; A.col = n; A.nnz = nnz;
62 
63  return A;
64 }
65 
81  const INT n,
82  const INT nnz)
83 {
84  iCSRmat A;
85 
86  if ( m > 0 ) {
87  A.IA = (INT *)fasp_mem_calloc(m+1, sizeof(INT));
88  }
89  else {
90  A.IA = NULL;
91  }
92 
93  if ( n > 0 ) {
94  A.JA = (INT *)fasp_mem_calloc(nnz, sizeof(INT));
95  }
96  else {
97  A.JA = NULL;
98  }
99 
100  if ( nnz > 0 ) {
101  A.val = (INT *)fasp_mem_calloc(nnz, sizeof(INT));
102  }
103  else {
104  A.val = NULL;
105  }
106 
107  A.row = m; A.col = n; A.nnz = nnz;
108 
109  return A;
110 }
111 
125 void fasp_dcsr_alloc (const INT m,
126  const INT n,
127  const INT nnz,
128  dCSRmat *A)
129 {
130  if ( m > 0 ) {
131  A->IA=(INT*)fasp_mem_calloc(m+1,sizeof(INT));
132  }
133  else {
134  A->IA = NULL;
135  }
136 
137  if ( n > 0 ) {
138  A->JA=(INT*)fasp_mem_calloc(nnz,sizeof(INT));
139  }
140  else {
141  A->JA = NULL;
142  }
143 
144  if ( nnz > 0 ) {
145  A->val=(REAL*)fasp_mem_calloc(nnz,sizeof(REAL));
146  }
147  else {
148  A->val = NULL;
149  }
150 
151  A->row=m; A->col=n; A->nnz=nnz;
152 
153  return;
154 }
155 
167 {
168  if ( A == NULL ) return;
169 
170  fasp_mem_free(A->IA); A->IA = NULL;
171  fasp_mem_free(A->JA); A->JA = NULL;
172  fasp_mem_free(A->val); A->val = NULL;
173 }
174 
186 {
187  if ( A == NULL ) return;
188 
189  fasp_mem_free(A->IA); A->IA = NULL;
190  fasp_mem_free(A->JA); A->JA = NULL;
191  fasp_mem_free(A->val); A->val = NULL;
192 }
193 
205 {
206  A->row = A->col = A->nnz = 0;
207  A->IA = A->JA = NULL;
208  A->val = NULL;
209 }
210 
222 {
223  A->row = A->col = A->nnz = 0;
224  A->IA = A->JA = NULL;
225  A->val = NULL;
226 }
227 
248  INT *P)
249 {
250  const INT n=A->row,nnz=A->nnz;
251  const INT *ia=A->IA, *ja=A->JA;
252  const REAL *Aval=A->val;
253  INT i,j,k,jaj,i1,i2,start;
254  INT nthreads = 1, use_openmp = FALSE;
255 
256 #ifdef _OPENMP
257  if ( MIN(n, nnz) > OPENMP_HOLDS ) {
258  use_openmp = TRUE;
259  nthreads = FASP_GET_NUM_THREADS();
260  }
261 #endif
262 
263  dCSRmat Aperm = fasp_dcsr_create(n,n,nnz);
264 
265  // form the transpose of P
266  INT *Pt = (INT*)fasp_mem_calloc(n,sizeof(INT));
267 
268  if (use_openmp) {
269  INT myid, mybegin, myend;
270 #ifdef _OPENMP
271 #pragma omp parallel for private(myid, mybegin, myend, i)
272 #endif
273  for (myid=0; myid<nthreads; ++myid) {
274  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
275  for (i=mybegin; i<myend; ++i) Pt[P[i]] = i;
276  }
277  }
278  else {
279  for (i=0; i<n; ++i) Pt[P[i]] = i;
280  }
281 
282  // compute IA of P*A (row permutation)
283  Aperm.IA[0] = 0;
284  for (i=0; i<n; ++i) {
285  k = P[i];
286  Aperm.IA[i+1] = Aperm.IA[i]+(ia[k+1]-ia[k]);
287  }
288 
289  // perform actual P*A
290  if (use_openmp) {
291  INT myid, mybegin, myend;
292 #ifdef _OPENMP
293 #pragma omp parallel for private(myid, mybegin, myend, i1, i2, k, start, j, jaj)
294 #endif
295  for (myid=0; myid<nthreads; ++myid) {
296  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
297  for (i=mybegin; i<myend; ++i) {
298  i1 = Aperm.IA[i]; i2 = Aperm.IA[i+1]-1;
299  k = P[i];
300  start = ia[k];
301  for (j=i1; j<=i2; ++j) {
302  jaj = start+j-i1;
303  Aperm.JA[j] = ja[jaj];
304  Aperm.val[j] = Aval[jaj];
305  }
306  }
307  }
308  }
309  else {
310  for (i=0; i<n; ++i) {
311  i1 = Aperm.IA[i]; i2 = Aperm.IA[i+1]-1;
312  k = P[i];
313  start = ia[k];
314  for (j=i1; j<=i2; ++j) {
315  jaj = start+j-i1;
316  Aperm.JA[j] = ja[jaj];
317  Aperm.val[j] = Aval[jaj];
318  }
319  }
320  }
321 
322  // perform P*A*P' (column permutation)
323  if (use_openmp) {
324  INT myid, mybegin, myend;
325 #ifdef _OPENMP
326 #pragma omp parallel for private(myid, mybegin, myend, k, j)
327 #endif
328  for (myid=0; myid<nthreads; ++myid) {
329  FASP_GET_START_END(myid, nthreads, nnz, &mybegin, &myend);
330  for (k=mybegin; k<myend; ++k) {
331  j = Aperm.JA[k];
332  Aperm.JA[k] = Pt[j];
333  }
334  }
335  }
336  else {
337  for (k=0; k<nnz; ++k) {
338  j = Aperm.JA[k];
339  Aperm.JA[k] = Pt[j];
340  }
341  }
342 
343  fasp_mem_free(Pt);
344 
345  return(Aperm);
346 }
347 
359 {
360  const INT n=A->col;
361  INT i,j,start,row_length;
362 
363  // temp memory for sorting rows of A
364  INT *index, *ja;
365  REAL *a;
366 
367  index = (INT*)fasp_mem_calloc(n,sizeof(INT));
368  ja = (INT*)fasp_mem_calloc(n,sizeof(INT));
369  a = (REAL*)fasp_mem_calloc(n,sizeof(REAL));
370 
371  for (i=0;i<n;++i) {
372  start=A->IA[i];
373  row_length=A->IA[i+1]-start;
374 
375  for (j=0;j<row_length;++j) index[j]=j;
376 
377  fasp_aux_iQuickSortIndex(&(A->JA[start]), 0, row_length-1, index);
378 
379  for (j=0;j<row_length;++j) {
380  ja[j]=A->JA[start+index[j]];
381  a[j]=A->val[start+index[j]];
382  }
383 
384  for (j=0;j<row_length;++j) {
385  A->JA[start+j]=ja[j];
386  A->val[start+j]=a[j];
387  }
388  }
389 
390  // clean up memory
391  fasp_mem_free(index);
392  fasp_mem_free(ja);
393  fasp_mem_free(a);
394 }
395 
411  dCSRmat *A,
412  dvector *diag)
413 {
414  INT i,k,j,ibegin,iend;
415 
416  if ( n==0 || n>A->row || n>A->col ) n = MIN(A->row,A->col);
417 
418  INT nthreads = 1, use_openmp = FALSE;
419 
420 #ifdef _OPENMP
421  if ( n > OPENMP_HOLDS ) {
422  use_openmp = TRUE;
423  nthreads = FASP_GET_NUM_THREADS();
424  }
425 #endif
426 
427  fasp_dvec_alloc(n, diag);
428 
429  if (use_openmp) {
430  INT mybegin,myend,myid;
431 #ifdef _OPENMP
432 #pragma omp parallel for private(myid, mybegin, myend, i, ibegin, iend, k, j)
433 #endif
434  for (myid = 0; myid < nthreads; myid++ ) {
435  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
436  for (i=mybegin; i<myend; i++) {
437  ibegin=A->IA[i]; iend=A->IA[i+1];
438  for (k=ibegin;k<iend;++k) {
439  j=A->JA[k];
440  if ((j-i)==0) {
441  diag->val[i] = A->val[k]; break;
442  } // end if
443  } // end for k
444  } // end for i
445  }
446  }
447  else {
448  for (i=0;i<n;++i) {
449  ibegin=A->IA[i]; iend=A->IA[i+1];
450  for (k=ibegin;k<iend;++k) {
451  j=A->JA[k];
452  if ((j-i)==0) {
453  diag->val[i] = A->val[k]; break;
454  } // end if
455  } // end for k
456  } // end for i
457  }
458 }
459 
474 void fasp_dcsr_getcol (const INT n,
475  dCSRmat *A,
476  REAL *col)
477 {
478  INT i,j, row_begin, row_end;
479  INT nrow = A->row, ncol = A->col;
480  INT status = FASP_SUCCESS;
481 
482  INT nthreads=1, use_openmp=FALSE;
483 
484 #ifdef _OPENMP
485  if ( nrow > OPENMP_HOLDS ) {
486  use_openmp = TRUE;
487  nthreads = FASP_GET_NUM_THREADS();
488  }
489 #endif
490 
491  // check the column index n
492  if ( n < 0 || n >= ncol ) {
493  printf("### ERROR: Column index %d is illegal!\n", n);
494  status = ERROR_DUMMY_VAR;
495  goto FINISHED;
496  }
497 
498  // get the column
499  if (use_openmp) {
500  INT mybegin, myend, myid;
501 
502 #ifdef _OPENMP
503 #pragma omp parallel for private(myid, mybegin, myend, i, j, row_begin, row_end)
504 #endif
505  for (myid = 0; myid < nthreads; myid++ ) {
506  FASP_GET_START_END(myid, nthreads, nrow, &mybegin, &myend);
507  for (i=mybegin; i<myend; i++) {
508  col[i] = 0.0;
509  row_begin = A->IA[i]; row_end = A->IA[i+1];
510  for (j=row_begin; j<row_end; ++j) {
511  if (A->JA[j] == n) {
512  col[i] = A->val[j];
513  }
514  } // end for j
515  } // end for i
516  }
517  }
518  else {
519  for (i=0; i<nrow; ++i) {
520  // set the entry to zero
521  col[i] = 0.0;
522  row_begin = A->IA[i]; row_end = A->IA[i+1];
523  for (j=row_begin; j<row_end; ++j) {
524  if (A->JA[j] == n) {
525  col[i] = A->val[j];
526  }
527  } // end for j
528  } // end for i
529  }
530 
531 FINISHED:
532  fasp_chkerr(status,__FUNCTION__);
533 }
534 
554 {
555  const INT num_rowsA = A -> row;
556  REAL * A_data = A -> val;
557  INT * A_i = A -> IA;
558  INT * A_j = A -> JA;
559 
560  // Local variable
561  INT i, j;
562  INT tempi, row_size;
563  REAL tempd;
564 
565 #ifdef _OPENMP
566  // variables for OpenMP
567  INT myid, mybegin, myend, ibegin, iend;
568  INT nthreads = FASP_GET_NUM_THREADS();
569 #endif
570 
571 #if DEBUG_MODE > 0
572  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
573 #endif
574 
575 #ifdef _OPENMP
576  if (num_rowsA > OPENMP_HOLDS) {
577 #pragma omp parallel for private(myid,i,j,ibegin,iend,tempi,tempd, mybegin, myend)
578  for (myid = 0; myid < nthreads; myid++) {
579  FASP_GET_START_END(myid, nthreads, num_rowsA, &mybegin, &myend);
580  for (i = mybegin; i < myend; i ++) {
581  ibegin = A_i[i]; iend = A_i[i+1];
582  // check whether the first entry is already diagonal
583  if (A_j[ibegin] != i) {
584  for (j = ibegin+1 ; j < iend; j ++) {
585  if (A_j[j] == i) {
586 #if DEBUG_MODE > 2
587  printf("### DEBUG: Switch entry_%d with entry_0\n", j);
588 #endif
589  tempi = A_j[ibegin];
590  A_j[ibegin] = A_j[j];
591  A_j[j] = tempi;
592 
593  tempd = A_data[ibegin];
594  A_data[ibegin] = A_data[j];
595  A_data[j] = tempd;
596  break;
597  }
598  }
599  if (j == iend) {
600  printf("### ERROR: Diagonal entry %d is missing or zero!\n", i);
601  fasp_chkerr(ERROR_MISC, __FUNCTION__);
602  }
603  }
604  }
605  }
606  }
607  else {
608 #endif
609  for (i = 0; i < num_rowsA; i ++) {
610  row_size = A_i[i+1] - A_i[i];
611  // check whether the first entry is already diagonal
612  if (A_j[0] != i) {
613  for (j = 1; j < row_size; j ++) {
614  if (A_j[j] == i) {
615 #if DEBUG_MODE > 2
616  printf("### DEBUG: Switch entry_%d with entry_0\n", j);
617 #endif
618  tempi = A_j[0];
619  A_j[0] = A_j[j];
620  A_j[j] = tempi;
621 
622  tempd = A_data[0];
623  A_data[0] = A_data[j];
624  A_data[j] = tempd;
625 
626  break;
627  }
628  }
629  if (j == row_size) {
630  printf("### ERROR: Diagonal entry %d is missing or zero!\n", i);
631  fasp_chkerr(ERROR_MISC, __FUNCTION__);
632  }
633  }
634  A_j += row_size;
635  A_data += row_size;
636  }
637 #ifdef _OPENMP
638  }
639 #endif
640 
641 #if DEBUG_MODE > 0
642  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
643 #endif
644 }
645 
660  REAL value)
661 {
662  const INT m = A->row;
663  const INT * ia = A->IA, * ja = A->JA;
664  REAL * aj = A->val;
665 
666  // Local variables
667  INT i,j,k,begin_row,end_row;
668  SHORT status=ERROR_UNKNOWN;
669 
670  for (i=0;i<m;++i) {
671  begin_row=ia[i],end_row=ia[i+1];
672  for (k=begin_row;k<end_row;++k) {
673  j=ja[k];
674  if (i==j) {
675  if (aj[k] < 0.0) goto FINISHED;
676  else if (aj[k] < SMALLREAL) aj[k]=value;
677  }
678  } // end for k
679  } // end for i
680 
681  status = FASP_SUCCESS;
682 
683 FINISHED:
684  return status;
685 }
686 
699  iCSRmat *B)
700 {
701  B->row=A->row;
702  B->col=A->col;
703  B->nnz=A->nnz;
704 
705  fasp_iarray_cp (A->row+1, A->IA, B->IA);
706  fasp_iarray_cp (A->nnz, A->JA, B->JA);
707  fasp_iarray_cp (A->nnz, A->val, B->val);
708 }
709 
724  dCSRmat *B)
725 {
726  B->row=A->row;
727  B->col=A->col;
728  B->nnz=A->nnz;
729 
730  fasp_iarray_cp (A->row+1, A->IA, B->IA);
731  fasp_iarray_cp (A->nnz, A->JA, B->JA);
732  fasp_array_cp (A->nnz, A->val, B->val);
733 }
734 
751  iCSRmat *AT)
752 {
753  const INT n=A->row, m=A->col, nnz=A->nnz, m1=m-1;
754 
755  // Local variables
756  INT i,j,k,p;
757  INT ibegin, iend;
758 
759 #if DEBUG_MODE > 1
760  printf("### DEBUG: m=%d, n=%d, nnz=%d\n",m,n,nnz);
761 #endif
762 
763  AT->row=m; AT->col=n; AT->nnz=nnz;
764 
765  AT->IA=(INT*)fasp_mem_calloc(m+1,sizeof(INT));
766 
767  AT->JA=(INT*)fasp_mem_calloc(nnz,sizeof(INT));
768 
769  if (A->val) {
770  AT->val=(INT*)fasp_mem_calloc(nnz,sizeof(INT));
771  }
772  else {
773  AT->val=NULL;
774  }
775 
776  // first pass: find the Number of nonzeros in the first m-1 columns of A
777  // Note: these Numbers are stored in the array AT.IA from 1 to m-1
778  fasp_iarray_set(m+1, AT->IA, 0);
779 
780  for (j=0;j<nnz;++j) {
781  i=A->JA[j]; // column Number of A = row Number of A'
782  if (i<m1) AT->IA[i+2]++;
783  }
784 
785  for (i=2;i<=m;++i) AT->IA[i]+=AT->IA[i-1];
786 
787  // second pass: form A'
788  if (A->val != NULL) {
789  for (i=0;i<n;++i) {
790  ibegin=A->IA[i], iend=A->IA[i+1];
791  for (p=ibegin;p<iend;p++) {
792  j=A->JA[p]+1;
793  k=AT->IA[j];
794  AT->JA[k]=i;
795  AT->val[k]=A->val[p];
796  AT->IA[j]=k+1;
797  } // end for p
798  } // end for i
799  }
800  else {
801  for (i=0;i<n;++i) {
802  ibegin=A->IA[i], iend=A->IA[i+1];
803  for (p=ibegin;p<iend;p++) {
804  j=A->JA[p]+1;
805  k=AT->IA[j];
806  AT->JA[k]=i;
807  AT->IA[j]=k+1;
808  } // end for p
809  } // end for i
810  } // end if
811 }
812 
827  dCSRmat *AT)
828 {
829  const INT n=A->row, m=A->col, nnz=A->nnz;
830 
831  // Local variables
832  INT i,j,k,p;
833 
834  AT->row=m;
835  AT->col=n;
836  AT->nnz=nnz;
837 
838  AT->IA=(INT*)fasp_mem_calloc(m+1,sizeof(INT));
839 
840  AT->JA=(INT*)fasp_mem_calloc(nnz,sizeof(INT));
841 
842  if (A->val) {
843  AT->val=(REAL*)fasp_mem_calloc(nnz,sizeof(REAL));
844 
845  }
846  else { AT->val=NULL; }
847 
848  // first pass: find the Number of nonzeros in the first m-1 columns of A
849  // Note: these Numbers are stored in the array AT.IA from 1 to m-1
850 
851  // fasp_iarray_set(m+1, AT->IA, 0);
852  memset(AT->IA, 0, sizeof(INT)*(m+1));
853 
854  for (j=0;j<nnz;++j) {
855  i=A->JA[j]; // column Number of A = row Number of A'
856  if (i<m-1) AT->IA[i+2]++;
857  }
858 
859  for (i=2;i<=m;++i) AT->IA[i]+=AT->IA[i-1];
860 
861  // second pass: form A'
862  if (A->val) {
863  for (i=0;i<n;++i) {
864  INT ibegin=A->IA[i], iend=A->IA[i+1];
865  for (p=ibegin;p<iend;p++) {
866  j=A->JA[p]+1;
867  k=AT->IA[j];
868  AT->JA[k]=i;
869  AT->val[k]=A->val[p];
870  AT->IA[j]=k+1;
871  } // end for p
872  } // end for i
873  }
874  else {
875  for (i=0;i<n;++i) {
876  INT ibegin=A->IA[i], iend1=A->IA[i+1];
877  for (p=ibegin;p<iend1;p++) {
878  j=A->JA[p]+1;
879  k=AT->IA[j];
880  AT->JA[k]=i;
881  AT->IA[j]=k+1;
882  } // end for p
883  } // end of i
884  } // end if
885 
886  return FASP_SUCCESS;
887 }
888 
889 /*
890  * \fn void fasp_dcsr_transpose (INT *row[2], INT *col[2], REAL *val[2],
891  * INT *nn, INT *tniz)
892  *
893  * \brief Transpose of a CSR matrix
894  *
895  * \param row[2] Pointers of the rows of the matrix and its transpose
896  * \param col[2] Pointers of the columns of the matrix and its transpose
897  * \param val[2] Pointers to the values of the matrix and its transpose
898  * \param nn Number of rows and columns of the matrices A and A'
899  * \param tniz Number of the nonzeros in the matrices A and A'
900  *
901  * \author Shuo Zhang
902  * \date 07/06/2009
903  *
904  * \note This subroutine transpose in CSR format IN ORDER
905  */
906 void fasp_dcsr_transpose (INT *row[2],
907  INT *col[2],
908  REAL *val[2],
909  INT *nn,
910  INT *tniz)
911 {
912  const INT nca=nn[1]; // Number of columns
913 
914  INT *izc = (INT *)fasp_mem_calloc(nn[1],sizeof(INT));
915  INT *izcaux = (INT *)fasp_mem_calloc(nn[1],sizeof(INT));
916 
917  // Local variables
918  INT i,m,itmp;
919 
920  // first pass: to set order right
921  for (i=0;i<tniz[0];++i) izc[col[0][i]]++;
922 
923  izcaux[0]=0;
924  for (i=1;i<nca;++i) izcaux[i]=izcaux[i-1]+izc[i-1];
925 
926  // second pass: form transpose
927  memset(izc,0,nca*sizeof(INT));
928 
929  for (i=0;i<tniz[0];++i) {
930  m=col[0][i];
931  itmp=izcaux[m]+izc[m];
932  row[1][itmp]=m;
933  col[1][itmp]=row[0][i];
934  val[1][itmp]=val[0][i];
935  izc[m]++;
936  }
937 
938  fasp_mem_free(izc);
939  fasp_mem_free(izcaux);
940 }
941 
958  dCSRmat *B,
959  REAL dtol)
960 {
961  INT i, j, k;
962  INT ibegin,iend1;
963 
964  INT nthreads = 1, use_openmp = FALSE;
965 
966 #ifdef _OPENMP
967  if ( B->nnz > OPENMP_HOLDS) {
968  use_openmp = TRUE;
969  nthreads = FASP_GET_NUM_THREADS();
970  }
971 #endif
972 
973  INT *index=(INT*)fasp_mem_calloc(A->nnz,sizeof(INT));
974 
975  B->row=A->row; B->col=A->col;
976 
977  B->IA=(INT*)fasp_mem_calloc(A->row+1,sizeof(INT));
978 
979  B->IA[0]=A->IA[0];
980 
981  // first pass: determine the size of B
982  k=0;
983  for (i=0;i<A->row;++i) {
984  ibegin=A->IA[i]; iend1=A->IA[i+1];
985  for (j=ibegin;j<iend1;++j)
986  if (ABS(A->val[j])>dtol) {
987  index[k]=j;
988  ++k;
989  } /* end of j */
990  B->IA[i+1]=k;
991  } /* end of i */
992  B->nnz=k;
993 
994  B->JA=(INT*)fasp_mem_calloc(B->nnz,sizeof(INT));
995  B->val=(REAL*)fasp_mem_calloc(B->nnz,sizeof(REAL));
996 
997  // second pass: generate the index and element to B
998  if ( use_openmp ) {
999  INT myid, mybegin, myend;
1000 #ifdef _OPENMP
1001 #pragma omp parallel for private(myid, i, mybegin, myend)
1002 #endif
1003  for (myid=0; myid<nthreads; myid++) {
1004  FASP_GET_START_END(myid, nthreads, B->nnz, &mybegin, &myend);
1005  for (i=mybegin; i<myend; ++i) {
1006  B->JA[i]=A->JA[index[i]];
1007  B->val[i]=A->val[index[i]];
1008  }
1009  }
1010  }
1011  else {
1012  for (i=0;i<B->nnz;++i) {
1013  B->JA[i]=A->JA[index[i]];
1014  B->val[i]=A->val[index[i]];
1015  }
1016  }
1017 
1018  fasp_mem_free(index);
1019 }
1020 
1038  REAL dtol)
1039 {
1040  const INT row=A->row;
1041  const INT nnz=A->nnz;
1042 
1043  INT i, j, k;
1044  INT ibegin, iend = A->IA[0];
1045  SHORT status = FASP_SUCCESS;
1046 
1047  k = 0;
1048  for ( i=0; i<row; ++i ) {
1049  ibegin = iend; iend = A->IA[i+1];
1050  for ( j=ibegin; j<iend; ++j )
1051  if ( ABS(A->val[j]) > dtol ) {
1052  A->JA[k] = A->JA[j];
1053  A->val[k] = A->val[j];
1054  ++k;
1055  } /* end of j */
1056  A->IA[i+1] = k;
1057  } /* end of i */
1058 
1059  if ( k <= nnz ) {
1060  A->nnz=k;
1061  A->JA = (INT *)fasp_mem_realloc(A->JA, k*sizeof(INT));
1062  A->val = (REAL *)fasp_mem_realloc(A->val, k*sizeof(REAL));
1063  }
1064  else {
1065  printf("### ERROR: Size of compressed matrix is larger than the original!\n");
1066  status = ERROR_UNKNOWN;
1067  }
1068 
1069  return (status);
1070 }
1071 
1086  INT offset)
1087 {
1088  const INT nnz=A->nnz;
1089  const INT n=A->row+1;
1090  INT i, *ai=A->IA, *aj=A->JA;
1091  INT nthreads = 1, use_openmp = FALSE;
1092 
1093 #ifdef _OPENMP
1094  if ( MIN(n, nnz) > OPENMP_HOLDS ) {
1095  use_openmp = TRUE;
1096  nthreads = FASP_GET_NUM_THREADS();
1097  }
1098 #endif
1099 
1100  if (use_openmp) {
1101  INT myid, mybegin, myend;
1102 #ifdef _OPENMP
1103 #pragma omp parallel for private(myid, mybegin, myend, i)
1104 #endif
1105  for (myid=0; myid<nthreads; myid++) {
1106  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
1107  for (i=mybegin; i<myend; i++) {
1108  ai[i] += offset;
1109  }
1110  }
1111  }
1112  else {
1113  for (i=0; i<n; ++i) ai[i]+=offset;
1114  }
1115 
1116  if (use_openmp) {
1117  INT myid, mybegin, myend;
1118 #ifdef _OPENMP
1119 #pragma omp parallel for private(myid, mybegin, myend, i)
1120 #endif
1121  for (myid=0; myid<nthreads; myid++) {
1122  FASP_GET_START_END(myid, nthreads, nnz, &mybegin, &myend);
1123  for (i=mybegin; i<myend; i++) {
1124  aj[i] += offset;
1125  }
1126  }
1127  }
1128  else {
1129  for (i=0; i<nnz; ++i) aj[i]+=offset;
1130  }
1131 }
1132 
1147  dvector *diag)
1148 {
1149  // information about matrix A
1150  const INT n = A->row;
1151  INT *IA = A->IA;
1152  INT *JA = A->JA;
1153  REAL *val = A->val;
1154 
1155  INT nthreads = 1, use_openmp = FALSE;
1156 
1157 #ifdef _OPENMP
1158  if ( n > OPENMP_HOLDS ) {
1159  use_openmp = TRUE;
1160  nthreads = FASP_GET_NUM_THREADS();
1161  }
1162 #endif
1163 
1164  // local variables
1165  INT i, j, k, row_start, row_end;
1166 
1167  if (diag->row != n) {
1168  printf("### ERROR: Size of diag = %d != size of matrix = %d!", diag->row, n);
1169  fasp_chkerr(ERROR_MISC, __FUNCTION__);
1170  }
1171 
1172  // work space
1173  REAL *work = (REAL *)fasp_mem_calloc(n, sizeof(REAL));
1174 
1175  if (use_openmp) {
1176  INT myid, mybegin, myend;
1177 #ifdef _OPENMP
1178 #pragma omp parallel for private(myid, mybegin, myend, i)
1179 #endif
1180  for (myid=0; myid<nthreads; myid++) {
1181  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
1182  for (i=mybegin; i<myend; i++) work[i] = sqrt(diag->val[i]);
1183  }
1184  }
1185  else {
1186  // square root of diagonal entries
1187  for (i=0; i<n; i++) work[i] = sqrt(diag->val[i]);
1188  }
1189 
1190  if (use_openmp) {
1191  INT myid, mybegin, myend;
1192 #ifdef _OPENMP
1193 #pragma omp parallel for private(myid, mybegin, myend, row_start, row_end, i, j, k)
1194 #endif
1195  for (myid=0; myid<nthreads; myid++) {
1196  FASP_GET_START_END(myid, nthreads, n, &mybegin, &myend);
1197  for (i=mybegin; i<myend; i++) {
1198  row_start = IA[i]; row_end = IA[i+1];
1199  for (j=row_start; j<row_end; j++) {
1200  k = JA[j];
1201  val[j] = val[j]/(work[i]*work[k]);
1202  }
1203  }
1204  }
1205  }
1206  else {
1207  // main loop
1208  for (i=0; i<n; i++) {
1209  row_start = IA[i]; row_end = IA[i+1];
1210  for (j=row_start; j<row_end; j++) {
1211  k = JA[j];
1212  val[j] = val[j]/(work[i]*work[k]);
1213  }
1214  }
1215  }
1216 
1217  // free work space
1218  if (work) fasp_mem_free(work);
1219 }
1220 
1233 {
1234  //local variable
1235  dCSRmat AT;
1236 
1237  //return variable
1238  dCSRmat SA;
1239 
1240  // get the transpose of A
1241  fasp_dcsr_trans (A, &AT);
1242 
1243  // get symmetrized A
1244  fasp_blas_dcsr_add (A, 1.0, &AT, 0.0, &SA);
1245 
1246  // clean
1247  fasp_dcsr_free(&AT);
1248 
1249  // return
1250  return SA;
1251 }
1252 
1266  INT *flags,
1267  INT *groups)
1268 {
1269 #ifdef MULTI_COLOR_ORDER
1270  INT k,i,j,pre,group;
1271  INT iend;
1272  INT icount;
1273  INT front,rear;
1274  INT n=A->row;
1275  INT *IA = A -> IA;
1276  INT *JA = A -> JA;
1277  INT *cq = (INT *)malloc(sizeof(INT)*(n+1));
1278  INT *newr = (INT *)malloc(sizeof(INT)*(n+1));
1279 
1280 #ifdef _OPENMP
1281 #pragma omp parallel for private(k)
1282 #endif
1283  for (k=0;k<n;k++) cq[k]= k;
1284 
1285  group = 0;
1286  for (k=0;k<n;k++) {
1287  if ((IA[k+1]-IA[k]) > group )
1288  group = IA[k+1]-IA[k];
1289  }
1290 
1291  A->IC = (INT *)malloc(sizeof(INT)*(group+2));
1292  A->ICMAP = (INT *)malloc(sizeof(INT)*(n));
1293 
1294  front = n-1; rear = n-1;
1295 
1296  memset(newr, -1, sizeof(INT)*(n+1));
1297  memset(A->ICMAP, 0, sizeof(INT)*n);
1298 
1299  group = 0; icount = 0;
1300  A->IC[ 0 ] = 0; pre=0;
1301 
1302  do {
1303  front ++; if (front == n) front =0;
1304  i = cq[front];
1305  if (i <= pre) {
1306  A->IC[group] = icount; A->ICMAP[icount] = i;
1307  group++; icount++;
1308  iend = IA[i+1];
1309  for (j= IA[i]; j< iend; j++) newr[ JA[j] ] = group;
1310  }
1311  else if (newr[i] == group) {
1312  rear ++; if (rear == n) rear = 0;
1313  cq[rear] = i;
1314  }
1315  else {
1316  A->ICMAP[icount] = i;
1317  icount++;
1318  iend = IA[i+1];
1319  for (j = IA[i]; j< iend; j++) newr[ JA[j] ] = group;
1320  }
1321  pre=i;
1322 
1323  } while (rear != front);
1324 
1325  A->IC[group] = icount;
1326  A->color = group;
1327  free(cq);
1328  free(newr);
1329  *groups = group;
1330 #else
1331  printf("### ERROR: MULTI_COLOR_ORDER has not been defined!\n");
1332 #endif
1333 }
1334 
1367  INT *p,
1368  dCSRmat *AT)
1369 {
1370  /* tested for permutation and transposition */
1371  /* transpose or permute; if A.val is null ===> transpose the
1372  structure only */
1373  const INT n=A->row, m=A->col, nnz=A->nnz;
1374  const INT *ia=NULL,*ja=NULL;
1375  const REAL *a=NULL;
1376  INT m1=m+1;
1377  ia=A->IA; ja=A->JA; a=A->val;
1378  /* introducing few extra pointers hould not hurt too much the speed */
1379  INT *iat=NULL, *jat=NULL;
1380  REAL *at=NULL;
1381 
1382  /* loop variables */
1383  INT i,j,jp,pi,iabeg,iaend,k;
1384 
1385  /* initialize */
1386  AT->row=m; AT->col=n; AT->nnz=nnz;
1387 
1388  /* all these should be allocated or change this to allocate them here */
1389  iat=AT->IA; jat=AT->JA; at=AT->val;
1390  for (i = 0; i < m1; ++i) iat[i] = 0;
1391  iaend=ia[n];
1392  for (i = 0; i < iaend; ++i) {
1393  j = ja[i] + 2;
1394  if (j < m1) iat[j]++;
1395  }
1396  iat[0] = 0;
1397  iat[1] = 0;
1398  if (m != 1) {
1399  for (i= 2; i < m1; ++i) {
1400  iat[i] += iat[i-1];
1401  }
1402  }
1403 
1404  if (p && a) {
1405  /* so we permute and also use matrix entries */
1406  for (i = 0; i < n; ++i) {
1407  pi=p[i];
1408  iabeg = ia[pi];
1409  iaend = ia[pi+1];
1410  if (iaend > iabeg){
1411  for (jp = iabeg; jp < iaend; ++jp) {
1412  j = ja[jp]+1;
1413  k = iat[j];
1414  jat[k] = i;
1415  at[k] = a[jp];
1416  iat[j] = k+1;
1417  }
1418  }
1419  }
1420  } else if (a && !p) {
1421  /* transpose values, no permutation */
1422  for (i = 0; i < n; ++i) {
1423  iabeg = ia[i];
1424  iaend = ia[i+1];
1425  if (iaend > iabeg){
1426  for (jp = iabeg; jp < iaend; ++jp) {
1427  j = ja[jp]+1;
1428  k = iat[j];
1429  jat[k] = i;
1430  at[k] = a[jp];
1431  iat[j] = k+1;
1432  }
1433  }
1434  }
1435  } else if (!a && p) {
1436  /* Only integers and permutation (only a is null) */
1437  for (i = 0; i < n; ++i) {
1438  pi=p[i];
1439  iabeg = ia[pi];
1440  iaend = ia[pi+1];
1441  if (iaend > iabeg){
1442  for (jp = iabeg; jp < iaend; ++jp) {
1443  j = ja[jp]+1;
1444  k = iat[j];
1445  jat[k] = i;
1446  iat[j] = k+1;
1447  }
1448  }
1449  }
1450  } else {
1451  /* Only integers and no permutation (both a and p are null */
1452  for (i = 0; i < n; ++i) {
1453  iabeg = ia[i];
1454  iaend = ia[i+1];
1455  if (iaend > iabeg){
1456  for (jp = iabeg; jp < iaend; ++jp) {
1457  j = ja[jp]+1;
1458  k = iat[j];
1459  jat[k] = i;
1460  iat[j] = k+1;
1461  }
1462  }
1463  }
1464  }
1465 
1466  return;
1467 }
1468 
1487  INT *p)
1488 {
1489  const INT n = A->row, nnz = A->nnz;
1490  dCSRmat Aperm1, Aperm;
1491 
1492  Aperm1 = fasp_dcsr_create(n,n,nnz);
1493  Aperm = fasp_dcsr_create(n,n,nnz);
1494 
1495  fasp_dcsr_transz(A,p,&Aperm1);
1496  fasp_dcsr_transz(&Aperm1,p,&Aperm);
1497 
1498  // clean up
1499  fasp_dcsr_free(&Aperm1);
1500 
1501  return(Aperm);
1502 }
1503 
1519  const SHORT isym)
1520 {
1521  const INT n = A->row, m = A->col, nnz = A->nnz;
1522  dCSRmat AT = fasp_dcsr_create(m,n,nnz);
1523 
1524  /* watch carefully who is a pointer and who is not in fasp_dcsr_transz() */
1525  fasp_dcsr_transz(A, NULL , &AT);
1526 
1527  /* if the matrix is symmetric, then only one transpose is needed
1528  and now we just copy */
1529  if ((m==n) && (isym))
1530  fasp_dcsr_cp(&AT, A);
1531  else
1532  fasp_dcsr_transz(&AT, NULL , A);
1533 
1534  // clean up
1535  fasp_dcsr_free(&AT);
1536 }
1537 
1538 /*---------------------------------*/
1539 /*-- End of File --*/
1540 /*---------------------------------*/
#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
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: message.c:199
void fasp_dcsr_diagpref(dCSRmat *A)
Re-order the column and data arrays of a CSR matrix, so that the first entry in each row is the diago...
Definition: sparse_csr.c:553
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
void fasp_dcsr_symdiagscale(dCSRmat *A, dvector *diag)
Symmetric diagonal scaling D^{-1/2}AD^{-1/2}.
Definition: sparse_csr.c:1146
#define REAL
Definition: fasp.h:67
void fasp_aux_iQuickSortIndex(INT *a, INT left, INT right, INT *index)
Reorder the index of (INT type) so that 'a' is in ascending order.
Definition: ordering.c:279
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 fasp_dcsr_trans(dCSRmat *A, dCSRmat *AT)
Find transpose of dCSRmat matrix A.
Definition: sparse_csr.c:826
void fasp_dcsr_null(dCSRmat *A)
Initialize CSR sparse matrix.
Definition: sparse_csr.c:204
#define ERROR_UNKNOWN
Definition: fasp_const.h:62
REAL * val
nonzero entries of A
Definition: fasp.h:166
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
Vector with n entries of REAL type.
Definition: fasp.h:342
void fasp_dcsr_multicoloring(dCSRmat *A, INT *flags, INT *groups)
Use the greedy multi-coloring to get color groups of the adjacency graph of A.
Definition: sparse_csr.c:1265
#define INT
Definition: fasp.h:64
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:193
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
INT row
row number of matrix A, m
Definition: fasp.h:181
void fasp_dcsr_sortz(dCSRmat *A, const SHORT isym)
Sort each row of A in ascending order w.r.t. column indices.
Definition: sparse_csr.c:1518
void FASP_GET_START_END(INT procid, INT nprocs, INT n, INT *start, INT *end)
Assign Load to each thread.
Definition: threads.c:83
INT fasp_blas_dcsr_add(dCSRmat *A, const REAL alpha, dCSRmat *B, const REAL beta, dCSRmat *C)
compute C = alpha*A + beta*B in CSR format
Definition: blas_csr.c:48
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_icsr_cp(iCSRmat *A, iCSRmat *B)
Copy a iCSRmat to a new one B=A.
Definition: sparse_csr.c:698
void fasp_iarray_cp(const INT n, INT *x, INT *y)
Copy an array to the other y=x.
Definition: array.c:185
void fasp_dcsr_cp(dCSRmat *A, dCSRmat *B)
copy a dCSRmat to a new one B=A
Definition: sparse_csr.c:723
void fasp_dcsr_transz(dCSRmat *A, INT *p, dCSRmat *AT)
Generalized transpose of A: (n x m) matrix given in dCSRmat format.
Definition: sparse_csr.c:1366
void fasp_dcsr_shift(dCSRmat *A, INT offset)
Re-index a REAL matrix in CSR format to make the index starting from 0 or 1.
Definition: sparse_csr.c:1085
#define OPENMP_HOLDS
Definition: fasp_const.h:248
INT nnz
number of nonzero entries
Definition: fasp.h:157
void fasp_icsr_null(iCSRmat *A)
Initialize CSR sparse matrix.
Definition: sparse_csr.c:221
void fasp_dcsr_getcol(const INT n, dCSRmat *A, REAL *col)
Get the n-th column of a CSR matrix A.
Definition: sparse_csr.c:474
void fasp_icsr_trans(iCSRmat *A, iCSRmat *AT)
Find transpose of iCSRmat matrix A.
Definition: sparse_csr.c:750
#define ERROR_MISC
Definition: fasp_const.h:35
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:190
void * fasp_mem_realloc(void *oldmem, LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: memory.c:110
INT * val
nonzero entries of A
Definition: fasp.h:196
iCSRmat fasp_icsr_create(const INT m, const INT n, const INT nnz)
Create CSR sparse matrix data memory space.
Definition: sparse_csr.c:80
SHORT fasp_dcsr_compress_inplace(dCSRmat *A, REAL dtol)
Compress a CSR matrix A IN PLACE by dropping small entries abs(aij)<=dtol.
Definition: sparse_csr.c:1037
void fasp_dcsr_compress(dCSRmat *A, dCSRmat *B, REAL dtol)
Compress a CSR matrix A and store in CSR matrix B by dropping small entries abs(aij)<=dtol.
Definition: sparse_csr.c:957
dCSRmat fasp_dcsr_permz(dCSRmat *A, INT *p)
Permute rows and cols of A, i.e. A=PAP' by the ordering in p.
Definition: sparse_csr.c:1486
INT row
row number of matrix A, m
Definition: fasp.h:151
void fasp_dcsr_sort(dCSRmat *A)
Sort each row of A in ascending order w.r.t. column indices.
Definition: sparse_csr.c:358
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
dCSRmat fasp_dcsr_perm(dCSRmat *A, INT *P)
Apply permutation of A, i.e. Aperm=PAP' by the orders given in P.
Definition: sparse_csr.c:247
void fasp_dcsr_getdiag(INT n, dCSRmat *A, dvector *diag)
Get first n diagonal entries of a CSR matrix A.
Definition: sparse_csr.c:410
#define ERROR_DUMMY_VAR
Definition: fasp_const.h:40
INT row
number of rows
Definition: fasp.h:345
void fasp_icsr_free(iCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: sparse_csr.c:185
#define ABS(a)
Definition: fasp.h:74
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define MIN(a, b)
Definition: fasp.h:73
void fasp_dvec_alloc(const INT m, dvector *u)
Create dvector data space of REAL type.
Definition: vec.c:99
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
INT nnz
number of nonzero entries
Definition: fasp.h:187
#define FALSE
Definition: fasp_const.h:68
INT col
column of matrix A, n
Definition: fasp.h:184
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
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
SHORT fasp_dcsr_regdiag(dCSRmat *A, REAL value)
Regularize diagonal entries of a CSR sparse matrix.
Definition: sparse_csr.c:659
Sparse matrix of INT type in CSR format.
Definition: fasp.h:178
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: sparse_csr.c:166
#define SMALLREAL
Definition: fasp_const.h:238
dCSRmat fasp_dcsr_sympat(dCSRmat *A)
Get symmetric part of a dCSRmat matrix.
Definition: sparse_csr.c:1232