Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
sparse_block.c
Go to the documentation of this file.
1 
6 #include <time.h>
7 
8 #ifdef _OPENMP
9 #include <omp.h>
10 #endif
11 
12 #include "fasp.h"
13 #include "fasp_block.h"
14 #include "fasp_functs.h"
15 
16 /*---------------------------------*/
17 /*-- Public Functions --*/
18 /*---------------------------------*/
19 
31 {
32  if (A == NULL) return; // Nothing need to be freed!
33 
34  INT i;
35  INT num_blocks = (A->brow)*(A->bcol);
36 
37  for ( i=0; i<num_blocks; i++ ) {
38  fasp_dcsr_free(A->blocks[i]);
39  A->blocks[i] = NULL;
40  }
41 
43  A->blocks = NULL;
44 }
45 
67  INT *Is,
68  INT *Js,
69  const INT m,
70  const INT n,
71  dCSRmat *B)
72 {
73  INT status = FASP_SUCCESS;
74 
75  INT i,j,k,nnz=0;
76  INT *col_flag;
77 
78  INT use_openmp = FALSE;
79 
80 #ifdef _OPENMP
81  INT stride_i, mybegin, myend, myid, nthreads;
82  if ( n > OPENMP_HOLDS ) {
83  use_openmp = TRUE;
84  nthreads = FASP_GET_NUM_THREADS();
85  }
86 #endif
87 
88  // create column flags
89  col_flag = (INT*)fasp_mem_calloc(A->col,sizeof(INT));
90 
91  B->row = m; B->col = n;
92 
93  B->IA = (INT*)fasp_mem_calloc(m+1,sizeof(INT));
94  B->JA = (INT*)fasp_mem_calloc(A->nnz,sizeof(INT));
95  B->val = (REAL*)fasp_mem_calloc(A->nnz,sizeof(REAL));
96 
97  if (use_openmp) {
98 #ifdef _OPENMP
99  stride_i = n/nthreads;
100 #pragma omp parallel private(myid, mybegin, myend, i) num_threads(nthreads)
101  {
102  myid = omp_get_thread_num();
103  mybegin = myid*stride_i;
104  if (myid < nthreads-1) myend = mybegin+stride_i;
105  else myend = n;
106  for (i=mybegin; i < myend; ++i) {
107  col_flag[Js[i]]=i+1;
108  }
109  }
110 #endif
111  }
112  else {
113  for (i=0;i<n;++i) col_flag[Js[i]]=i+1;
114  }
115 
116  // Count nonzeros for sub matrix and fill in
117  B->IA[0]=0;
118  for (i=0;i<m;++i) {
119  for (k=A->IA[Is[i]];k<A->IA[Is[i]+1];++k) {
120  j=A->JA[k];
121  if (col_flag[j]>0) {
122  B->JA[nnz]=col_flag[j]-1;
123  B->val[nnz]=A->val[k];
124  nnz++;
125  }
126  } /* end for k */
127  B->IA[i+1]=nnz;
128  } /* end for i */
129  B->nnz=nnz;
130 
131  // re-allocate memory space
132  B->JA=(INT*)fasp_mem_realloc(B->JA, sizeof(INT)*nnz);
133  B->val=(REAL*)fasp_mem_realloc(B->val, sizeof(REAL)*nnz);
134 
135  fasp_mem_free(col_flag);
136 
137  return(status);
138 }
139 
161  INT *Is,
162  INT *Js,
163  const INT m,
164  const INT n,
165  dBSRmat *B)
166 {
167  INT status = FASP_SUCCESS;
168  INT i,j,k,nnz=0;
169  INT *col_flag;
170  INT use_openmp = FALSE;
171 
172  const INT nb = A->nb;
173  const INT nb2=nb*nb;
174 
175 #ifdef _OPENMP
176  INT myid, mybegin, stride_i, myend, nthreads;
177  if ( n > OPENMP_HOLDS ) {
178  use_openmp = TRUE;
179  nthreads = FASP_GET_NUM_THREADS();
180  }
181 #endif
182 
183  // create colum flags
184  col_flag=(INT*)fasp_mem_calloc(A->COL,sizeof(INT));
185 
186  B->ROW=m; B->COL=n; B->nb=nb; B->storage_manner=A->storage_manner;
187 
188  B->IA=(INT*)fasp_mem_calloc(m+1,sizeof(INT));
189 
190  if (use_openmp) {
191 #ifdef _OPENMP
192  stride_i = n/nthreads;
193 #pragma omp parallel private(myid, mybegin, myend, i) num_threads(nthreads)
194  {
195  myid = omp_get_thread_num();
196  mybegin = myid*stride_i;
197  if (myid < nthreads-1) myend = mybegin+stride_i;
198  else myend = n;
199  for (i=mybegin; i < myend; ++i) {
200  col_flag[Js[i]]=i+1;
201  }
202  }
203 #endif
204  }
205  else {
206  for (i=0;i<n;++i) col_flag[Js[i]]=i+1;
207  }
208 
209  // first pass: count nonzeros for sub matrix
210  B->IA[0]=0;
211  for (i=0;i<m;++i) {
212  for (k=A->IA[Is[i]];k<A->IA[Is[i]+1];++k) {
213  j=A->JA[k];
214  if (col_flag[j]>0) nnz++;
215  } /* end for k */
216  B->IA[i+1]=nnz;
217  } /* end for i */
218  B->NNZ=nnz;
219 
220  // allocate
221  B->JA=(INT*)fasp_mem_calloc(nnz,sizeof(INT));
222 
223  B->val=(REAL*)fasp_mem_calloc(nnz*nb2,sizeof(REAL));
224 
225  // second pass: copy data to B
226  // TODO: No need to do the following loop, need to be modified!! --Xiaozhe
227  nnz = 0;
228  for (i=0;i<m;++i) {
229  for (k=A->IA[Is[i]];k<A->IA[Is[i]+1];++k) {
230  j=A->JA[k];
231  if (col_flag[j]>0) {
232  B->JA[nnz]=col_flag[j]-1;
233  memcpy(B->val+nnz*nb2, A->val+k*nb2, nb2*sizeof(REAL));
234  nnz++;
235  }
236  } /* end for k */
237  } /* end for i */
238 
239  fasp_mem_free(col_flag);
240 
241  return(status);
242 }
243 
257 {
258  // information about A
259  const INT ROW = A->ROW;
260  const INT COL = A->COL;
261  const INT NNZ = A->NNZ;
262  const SHORT nc = A->nb;
263  const INT nc2 = nc*nc;
264  const REAL TOL = 1e-8;
265 
266  REAL *val = A->val;
267  INT *IA = A->IA;
268  INT *JA = A->JA;
269 
270  // Pressure block
271  dCSRmat P_csr = fasp_dcsr_create(ROW, COL, NNZ);
272  REAL *Pval=P_csr.val;
273 
274  // get pressure block
275  memcpy(P_csr.JA, JA, NNZ*sizeof(INT));
276  memcpy(P_csr.IA, IA, (ROW+1)*sizeof(INT));
277 
278 #ifdef _OPENMP
279  INT i;
280 
281 #pragma omp parallel for if(NNZ>OPENMP_HOLDS)
282  for (i=NNZ-1; i>=0; i--) {
283  Pval[i] = val[i*nc2];
284  }
285 #else
286  INT i, j;
287 
288  for (i=NNZ, j=NNZ*nc2-nc2 + (0*nc+0); i--; j-=nc2) {
289  Pval[i] = val[j];
290  }
291 #endif
292 
293  // compress CSR format
294  fasp_dcsr_compress_inplace(&P_csr,TOL);
295 
296  // return P
297  return P_csr;
298 }
299 
313 {
314  // information about A
315  const INT ROW = A->ROW;
316  const INT COL = A->COL;
317  const INT NNZ = A->NNZ;
318  const SHORT nc = A->nb;
319  const INT nc2 = nc*nc;
320  const REAL TOL = 1e-8;
321 
322  REAL *val = A->val;
323  INT *IA = A->IA;
324  INT *JA = A->JA;
325 
326  // CSR matrix
327  dCSRmat Acsr = fasp_dcsr_create(ROW, COL, NNZ);
328  REAL *Aval=Acsr.val;
329 
330  // get structure
331  memcpy(Acsr.JA, JA, NNZ*sizeof(INT));
332  memcpy(Acsr.IA, IA, (ROW+1)*sizeof(INT));
333 
334  INT i, j, k;
335  INT row_start, row_end;
336 
337  for (i=0; i<=ROW; i++){
338 
339  row_start = A->IA[i]; row_end = A->IA[i+1];
340 
341  for (k = row_start; k<row_end; k++) {
342 
343  j = A->JA[k];
344 
345  if ( i == j ) {
346  Aval[k] = fasp_blas_smat_Linfinity(val+k*nc2, nc);
347  }
348  else
349  {
350  Aval[k] = (-1)*fasp_blas_smat_Linfinity(val+k*nc2, nc);
351  }
352 
353  }
354 
355  }
356 
357  // compress CSR format
358  fasp_dcsr_compress_inplace(&Acsr,TOL);
359 
360  // return CSR matrix
361  return Acsr;
362 }
363 
364 /*---------------------------------*/
365 /*-- End of File --*/
366 /*---------------------------------*/
#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
dCSRmat ** blocks
blocks of dCSRmat, point to blocks[brow][bcol]
Definition: fasp_block.h:93
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
#define REAL
Definition: fasp.h:67
INT nb
dimension of each sub-block
Definition: fasp_block.h:56
REAL * val
nonzero entries of A
Definition: fasp.h:166
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:44
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
#define INT
Definition: fasp.h:64
#define FASP_SUCCESS
Definition of return status and error messages.
Definition: fasp_const.h:27
INT bcol
column number of blocks A, n
Definition: fasp_block.h:90
dCSRmat fasp_dbsr_Linfinity_dcsr(dBSRmat *A)
get dCSRmat from a dBSRmat matrix using L_infinity norm of each small block
Definition: sparse_block.c:312
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
SHORT fasp_dbsr_getblk(dBSRmat *A, INT *Is, INT *Js, const INT m, const INT n, dBSRmat *B)
Get a sub BSR matrix of A with specified rows and columns.
Definition: sparse_block.c:160
#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.
void * fasp_mem_realloc(void *oldmem, LONGLONG tsize)
Reallocate, initiate, and check memory.
Definition: memory.c:110
SHORT fasp_dcsr_getblk(dCSRmat *A, INT *Is, INT *Js, const INT m, const INT n, dCSRmat *B)
Get a sub CSR matrix of A with specified rows and columns.
Definition: sparse_block.c:66
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
dCSRmat fasp_dbsr_getblk_dcsr(dBSRmat *A)
get dCSRmat block from a dBSRmat matrix
Definition: sparse_block.c:256
void fasp_bdcsr_free(block_dCSRmat *A)
Free block CSR sparse matrix data memory space.
Definition: sparse_block.c:30
INT row
row number of matrix A, m
Definition: fasp.h:151
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
INT brow
row number of blocks in A, m
Definition: fasp_block.h:87
REAL fasp_blas_smat_Linfinity(REAL *A, const INT n)
Compute the L infinity norm of A.
Definition: smat.c:595
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
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
#define FALSE
Definition: fasp_const.h:68
INT * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160
INT * IA
integer array of row pointers, the size is ROW+1
Definition: fasp_block.h:70
void fasp_dcsr_free(dCSRmat *A)
Free CSR sparse matrix data memory space.
Definition: sparse_csr.c:166