Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
rap.c
Go to the documentation of this file.
1 
11 #include <math.h>
12 #include <time.h>
13 
14 #include "fasp.h"
15 #include "fasp_functs.h"
16 
17 /*---------------------------------*/
18 /*-- Public Functions --*/
19 /*---------------------------------*/
20 
34  INT *jr,
35  REAL *r,
36  INT *ia,
37  INT *ja,
38  REAL *a,
39  INT *ipt,
40  INT *jpt,
41  REAL *pt,
42  INT n,
43  INT nc,
44  INT *maxrpout,
45  INT *ipin,
46  INT *jpin)
47 {
48  dCSRmat ac;
49  INT n1,nc1,nnzp,maxrp;
50  INT *ip=NULL,*jp=NULL;
51 
52  /*=========================================================*/
53  /*
54  if ipin is null, this
55  means that we need to do the transpose of p here; otherwise,
56  these are considered to be input
57  */
58  maxrp=0;
59  nnzp=ipt[nc]-1;
60  n1=n+1;
61 
62  if(!ipin) {
63  ip = (INT *)calloc(n1,sizeof(INT));
64  jp = (INT *)calloc(nnzp,sizeof(INT));
65  /* these must be null anyway, so no need to assign null
66  ipin=NULL;
67  jpin=NULL;
68  */
69  } else {
70  ip=ipin;
71  jp=jpin;
72  }
73 
74  fasp_sparse_iit_(ipt,jpt,&nc,&n,ip,jp);
75 
76  /* triple matrix product: R * A * transpose(P^T)=R*A*P.*/
77  /* A is square n by n*/
78  /* Note: to compute R*A* P the input are R, A and P^T */
79  /* we need to transpose now the structure of P, because the input is P^T */
80  /* end of transpose of the boolean corresponding to P */
81  /* ic are the addresses of the rows of the output */
82  nc1 = nc+1;
83  ac.IA = (INT *)calloc(nc1,sizeof(INT));
84 
85  /*
86  First call is with jc=null so that we find the number of
87  nonzeros in the result
88  */
89  ac.JA=NULL;
90  fasp_sparse_rapms_(ir,jr,ia,ja,ip,jp,&n,&nc,ac.IA,ac.JA,&maxrp);
91  ac.nnz = ac.IA[nc]-1;
92  ac.JA = (INT *)calloc(ac.nnz,sizeof(INT));
93 
94  /*
95  second call is to fill the column indexes array jc.
96  */
97  fasp_sparse_rapms_(ir,jr,ia,ja,ip,jp,&n,&nc,ac.IA,ac.JA,&maxrp);
98  if(!ipin){
99  if(ip) free(ip);
100  if(jp) free(jp);
101  }
102  ac.val = (REAL *)calloc(ac.nnz,sizeof(REAL));
103  /* this is the compute with the entries */
104  fasp_sparse_rapcmp_(ir,jr,r,ia,ja,a,ipt,jpt,pt,&n,&nc,ac.IA,ac.JA,ac.val,&maxrp);
105  ac.row=nc;
106  ac.col=nc;
107 
108  /*=========================================================*/
109  *maxrpout=maxrp;
110 
111  return ac;
112 }
113 
114 /*---------------------------------*/
115 /*-- End of File --*/
116 /*---------------------------------*/
INT * JA
integer array of column indexes, the size is nnz
Definition: fasp.h:163
#define REAL
Definition: fasp.h:67
void fasp_sparse_iit_(INT *ia, INT *ja, INT *na, INT *ma, INT *iat, INT *jat)
Transpose a boolean matrix (only given by ia, ja)
Definition: sparse_util.c:197
void fasp_sparse_rapms_(INT *ir, INT *jr, INT *ia, INT *ja, INT *ip, INT *jp, INT *nin, INT *ncin, INT *iac, INT *jac, INT *maxrout)
Calculates the nonzero structure of R*A*P, if jac is not null. If jac is null only finds num of nonze...
Definition: sparse_util.c:514
dCSRmat fasp_blas_dcsr_rap2(INT *ir, INT *jr, REAL *r, INT *ia, INT *ja, REAL *a, INT *ipt, INT *jpt, REAL *pt, INT n, INT nc, INT *maxrpout, INT *ipin, INT *jpin)
Compute R*A*P.
Definition: rap.c:33
REAL * val
nonzero entries of A
Definition: fasp.h:166
#define INT
Definition: fasp.h:64
INT nnz
number of nonzero entries
Definition: fasp.h:157
INT col
column of matrix A, n
Definition: fasp.h:154
Main header file for FASP.
void fasp_sparse_rapcmp_(INT *ir, INT *jr, REAL *r, INT *ia, INT *ja, REAL *a, INT *ipt, INT *jpt, REAL *pt, INT *nin, INT *ncin, INT *iac, INT *jac, REAL *ac, INT *idummy)
Calculates R*A*P after the nonzero structure of the result is known. iac,jac,ac have to be allocated ...
Definition: sparse_util.c:788
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 * IA
integer array of row pointers, the size is m+1
Definition: fasp.h:160