EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
evsl_cxsparse.c
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include "def.h"
5 #include "struct.h"
6 #include "CXSparse/Include/cs.h"
7 #include "internal_proto.h"
8 #include "evsl_direct.h"
9 
15 typedef struct _BSolDataDirect {
16  /* problem size */
17  int n;
18  /* symbolic factor */
20  /* numeric factor */
22  /* workspace for solve, of size n */
23  double *w;
25 
26 typedef struct _ASBSolDataDirect {
27  /* problem size */
28  int n;
29  /* symbolic factor */
31  /* numeric factor */
33  /* workspace for solve, of size n */
36 
37 
38 /* true for off-diagonal entries */
39 static int dropdiag_di (int i, int j, double aij, void *other) { return (i != j) ;}
40 
41 /* C = A + triu(A,1)' */
42 static cs_di *make_sym_di (cs_di *A)
43 {
44  cs_di *AT, *C ;
45  AT = cs_di_transpose (A, 1) ; /* AT = A' */
46  cs_di_fkeep (AT, &dropdiag_di, NULL) ; /* drop diagonal entries from AT */
47  C = cs_di_add (A, AT, 1, 1) ; /* C = A+AT */
48  cs_di_spfree (AT) ;
49  return (C) ;
50 }
51 
57 int SetupBSolDirect(csrMat *B, void **data) {
58  double tms = evsl_timer();
59  BSolDataDirect *Bsol_data;
60  Malloc(Bsol_data, 1, BSolDataDirect);
61  int i,j,k, n, nnz, *csp, *csi;
62  double *csx;
63  /* csparse csc matrix and the symmetrized one */
64  cs_di Bcs, *Bsym;
65  /* csparse factors: symbolic and numeric */
66  cs_dis *S;
67  cs_din *N;
68 
69  n = B->nrows;
70  nnz = B->ia[n];
71  /* allocate memory for Bcs */
72  Malloc(csp, n+1, int);
73  Malloc(csi, nnz, int);
74  Malloc(csx, nnz, double);
75  csp[0] = 0;
76  /* copy the lower part of B into cs */
77  /* nnz in cs */
78  k = 0;
79  for (i=0; i<n; i++) {
80  for (j=B->ia[i]; j<B->ia[i+1]; j++) {
81  int c = B->ja[j];
82  if (i >= c) {
83  csi[k] = c;
84  csx[k] = B->a[j];
85  k++;
86  }
87  }
88  csp[i+1] = k;
89  }
90  /* can view this as csc of upper part of B */
91  Bcs.nzmax = k;
92  Bcs.m = n;
93  Bcs.n = n;
94  Bcs.p = csp;
95  Bcs.i = csi;
96  Bcs.x = csx;
97  Bcs.nz = -1;
98  /* symmetrize B */
99  Bsym = make_sym_di(&Bcs);
100  if (!Bsym) {
101  return -1;
102  }
103  /* factorization of B */
104  S = cs_di_schol(1, Bsym) ; /* ordering and symbolic analysis */
105  N = cs_di_chol(Bsym, S) ; /* numeric Cholesky factorization */
106 
107  if (!(S && N)) {
108  return -2;
109  }
110 
111  //printf("Lnz %d %d\n", N->L->nzmax, N->L->p[N->L->n]);
112 
113  /* free */
114  cs_di_spfree(Bsym);
115  free(csp);
116  free(csi);
117  free(csx);
118 
119  Bsol_data->n = n;
120  Bsol_data->S = S;
121  Bsol_data->N = N;
122  Malloc(Bsol_data->w, n, double);
123 
124  *data = (void *) Bsol_data;
125 
126  double tme = evsl_timer();
127  evslstat.t_setBsv += tme - tms;
128 
129  return 0;
130 }
131 
135 void BSolDirect(double *b, double *x, void *data) {
136  BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
137  cs_dis *S = Bsol_data->S;
138  cs_din *N = Bsol_data->N;
139  int n = Bsol_data->n;
140  double *w = Bsol_data->w;
141 
142  cs_di_ipvec (S->pinv, b, w, n) ; /* w = P*b */
143  cs_di_lsolve (N->L, w) ; /* w = L\w */
144  cs_di_ltsolve (N->L, w) ; /* w = L'\w */
145  cs_di_pvec (S->pinv, w, x, n) ; /* x = P'*w */
146 }
147 
151 void LTSolDirect(double *b, double *x, void *data) {
152  BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
153  cs_dis *S = Bsol_data->S;
154  cs_din *N = Bsol_data->N;
155  int n = Bsol_data->n;
156  double *w = Bsol_data->w;
157 
158  memcpy(w, b, n*sizeof(double)); /* w = b */
159  cs_di_ltsolve (N->L, w) ; /* w = L'\w */
160  cs_di_pvec (S->pinv, w, x, n) ; /* x = P'*w */
161 }
162 
165 void FreeBSolDirectData(void *data) {
166  BSolDataDirect *Bsol_data = (BSolDataDirect *) data;
167  cs_dis *S = Bsol_data->S;
168  cs_din *N = Bsol_data->N;
169  double *w = Bsol_data->w;
170 
171  cs_di_sfree(S);
172  cs_di_nfree(N);
173  free(w);
174  free(data);
175 }
176 
190 int SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num,
191  complex double *zk, void **data) {
192  double tms = evsl_timer();
193  int i, j, nrow, ncol, nnzB, nnzC, *map;
194  csrMat *B, C, eye;
195  /* the shifted matrix
196  * C = A - s * B */
197  int *Cp, *Ci;
198  double *Cx, zkr1;
199  /* complex values of C */
200  cs_complex_t *Cz;
201  /* cs matrix complex integer */
202  cs_ci Ccs;
203  /* csparse factors: symbolic and numeric */
204  cs_cis *S = NULL;
205  cs_cin *N = NULL;
206  ASBSolDataDirect *ASBdata;
207 
208  nrow = A->nrows;
209  ncol = A->ncols;
210 
211  if (BB) {
212  B = BB;
213  } else {
214  /* if B==NULL, B=I, standard e.v. prob */
215  speye(nrow, &eye);
216  B = &eye;
217  }
218 
219  nnzB = B->ia[nrow];
220  /* map from nnz in B to nnz in C, useful for multi-poles */
221  Malloc(map, nnzB, int);
222  /* NOTE: if the matrix entries need to be sorted.
223  * The matadd routine will guarantee this
224  * C = A + 0.0 * B
225  * This actually computes the pattern of A + B
226  * and also can guarantee C has sorted rows
227  * map is the mapping from nnzB to nnzC */
228  matadd(1.0, 0.0, A, B, &C, NULL, map);
229  /* arrays for C */
230  nnzC = C.ia[nrow];
231  Cp = C.ia;
232  Ci = C.ja;
233  Cx = C.a;
234  /* complex values of C */
235  Malloc(Cz, nnzC, cs_complex_t);
236  for (i=0; i<nnzC; i++) {
237  Cz[i] = Cx[i] + 0.0 * I;
238  }
239  /* since C is complex symmetric, so its CSR and CSC are the same
240  * put them in a cs matrix */
241  Ccs.nzmax = nnzC;
242  Ccs.m = nrow;
243  Ccs.n = ncol;
244  Ccs.p = Cp;
245  Ccs.i = Ci;
246  Ccs.x = Cz;
247  Ccs.nz = -1;
248  /* pole loop
249  * for each pole we shift with B and factorize */
250  zkr1 = 0.0;
251  for (i=0; i<num; i++) {
252  /* the complex shift for pole i */
253  double zkr = creal(zk[i]);
254  double zkc = cimag(zk[i]);
255 
256  // shift B
257  for (j=0; j<nnzB; j++) {
258  int p = map[j];
259  double v = B->a[j];
260  CHKERR(Ci[p] != B->ja[j]);
261  double r = creal(Cz[p]) - (zkr - zkr1) * v;
262  Cz[p] = r - zkc * v * I;
263  }
264 
265  /* only do symbolic factorization once */
266  if (i == 0) {
267  /* Symbolic Factorization */
268  /* order == 1: sym, 0: not QR */
269  S = cs_ci_sqr(1, &Ccs, 0); /* ordering and symbolic analysis */
270  if (!S) {
271  return -1;
272  }
273  }
274  /* Numerical Factorization */
275  /* 1.0 results in conventional partial pivoting. */
276  /* 0.0 results in no pivoting. */
277  N = cs_ci_lu(&Ccs, S, 0.0);
278  if (!N) {
279  return -2;
280  }
281 
282  //printf("Lnz %d %d\n", N->L->nzmax, N->L->p[N->L->n]);
283  //printf("Unz %d %d\n", N->U->nzmax, N->U->p[N->U->n]);
284  //printf("%e\n", (N->L->nzmax + N->U->nzmax + 0.0) / nnzC);
285 
286  /* save the data */
287  Malloc(ASBdata, 1, ASBSolDataDirect);
288  ASBdata->n = nrow;
289  ASBdata->S = S;
290  ASBdata->N = N;
291  Malloc(ASBdata->b, nrow, cs_complex_t);
292  Malloc(ASBdata->x, nrow, cs_complex_t);
293  data[i] = ASBdata;
294 
295  /* for the next shift */
296  zkr1 = zkr;
297  } /* for (i=...)*/
298 
299  free(map);
300  free(Cz);
301  free_csr(&C);
302  if (!BB) {
303  free_csr(&eye);
304  }
305 
306  double tme = evsl_timer();
307  evslstat.t_setASigBsv += tme - tms;
308 
309  return 0;
310 }
311 
324 void ASIGMABSolDirect(int n, double *br, double *bi, double *xr,
325  double *xz, void *data) {
326 
327  ASBSolDataDirect *sol_data = (ASBSolDataDirect *) data;
328  int i;
329  CHKERR(n != sol_data->n);
330  cs_cis *S = sol_data->S;
331  cs_cin *N = sol_data->N;
332  cs_complex_t *b = sol_data->b;
333  cs_complex_t *x = sol_data->x;
334  /* copy rhs */
335  for (i=0; i<n; i++) {
336  b[i] = br[i] + bi[i] * I;
337  }
338  cs_ci_ipvec (N->pinv, b, x, n) ; /* x = b(p) */
339  cs_ci_lsolve (N->L, x) ; /* x = L\x */
340  cs_ci_usolve (N->U, x) ; /* x = U\x */
341  cs_ci_ipvec (S->q, x, b, n) ; /* b(q) = x */
342  /* copy sol */
343  for (i=0; i<n; i++) {
344  xr[i] = creal(b[i]);
345  xz[i] = cimag(b[i]);
346  }
347 }
348 
352 void FreeASIGMABSolDirect(int num, void **data) {
353  int i;
354  for (i=0; i<num; i++) {
355  ASBSolDataDirect *soldata = (ASBSolDataDirect *) data[i];
356  if (i == 0) {
357  cs_ci_sfree(soldata->S);
358  }
359  cs_ci_nfree(soldata->N);
360  free(soldata->b);
361  free(soldata->x);
362  free(soldata);
363  }
364 }
365 
int cs_ci_lsolve(const cs_ci *L, cs_complex_t *x)
void ASIGMABSolDirect(int n, double *br, double *bi, double *xr, double *xz, void *data)
complex linear solver routine passed to evsl
int m
Definition: cs.h:66
int * p
Definition: cs.h:342
defs in EVSL
int SetupASIGMABSolDirect(csrMat *A, csrMat *BB, int num, complex double *zk, void **data)
setup CXsparse solver for A - SIGMA B
cs_di * cs_di_add(const cs_di *A, const cs_di *B, double alpha, double beta)
int * q
Definition: cs.h:377
int ncols
Definition: struct.h:32
#define CHKERR(ierr)
Definition: def.h:19
int SetupBSolDirect(csrMat *B, void **data)
Setup the B-sol by computing the Cholesky factorization of B.
Definition: evsl_cxsparse.c:57
cs_complex_t * x
Definition: cs.h:344
void FreeASIGMABSolDirect(int num, void **data)
free the data needed by CXSparse
int * pinv
Definition: cs.h:101
cs_di * cs_di_transpose(const cs_di *A, int values)
int m
Definition: cs.h:340
double * x
Definition: cs.h:70
#define cs_complex_t
Definition: cs.h:41
int cs_di_pvec(const int *p, const double *b, double *x, int n)
int cs_di_fkeep(cs_di *A, int(*fkeep)(int, int, double, void *), void *other)
void FreeBSolDirectData(void *data)
Free solver data.
void free_csr(csrMat *csr)
memory deallocation for csr matrix
Definition: spmat.c:91
This file contains function prototypes and constant definitions internally used in EVSL...
cs_cin * cs_ci_lu(const cs_ci *A, const cs_cis *S, double tol)
int n
Definition: cs.h:67
int nz
Definition: cs.h:345
int * pinv
Definition: cs.h:390
int cs_ci_ipvec(const int *p, const cs_complex_t *b, cs_complex_t *x, int n)
int speye(int n, csrMat *A)
return an identity matrix of dimension n
Definition: spmat.c:357
#define Malloc(base, nmem, type)
Definition: def.h:22
int cs_ci_usolve(const cs_ci *U, cs_complex_t *x)
cs_din * cs_di_nfree(cs_din *N)
evslStat evslstat
global statistics of EVSL
Definition: evsl.c:21
int * i
Definition: cs.h:69
struct _ASBSolDataDirect ASBSolDataDirect
structs used in evsl
cs_ci * L
Definition: cs.h:388
int * i
Definition: cs.h:343
sparse matrix format: the compressed sparse row (CSR) format, 0-based
Definition: struct.h:31
double * a
Definition: struct.h:37
int nz
Definition: cs.h:71
int nzmax
Definition: cs.h:339
int cs_di_ipvec(const int *p, const double *b, double *x, int n)
int * ia
Definition: struct.h:32
void LTSolDirect(double *b, double *x, void *data)
Solver function of L^{T} x = L^{-T}*b.
int n
Definition: cs.h:341
double t_setBsv
Definition: struct.h:198
int matadd(double alp, double bet, csrMat *A, csrMat *B, csrMat *C, int *mapA, int *mapB)
matrix addition C = alp * A + bet * B
Definition: spmat.c:307
cs_dis * cs_di_schol(int order, const cs_di *A)
int nzmax
Definition: cs.h:65
int * ja
Definition: struct.h:32
cs_din * cs_di_chol(const cs_di *A, const cs_dis *S)
int cs_di_ltsolve(const cs_di *L, double *x)
struct _BSolDataDirect BSolDataDirect
cs_complex_t * b
Definition: evsl_cxsparse.c:34
cs_ci * U
Definition: cs.h:389
cs_cis * cs_ci_sfree(cs_cis *S)
double t_setASigBsv
Definition: struct.h:199
void BSolDirect(double *b, double *x, void *data)
Solver function of B.
Definitions used for direct solver interface.
cs_di * L
Definition: cs.h:113
int nrows
Definition: struct.h:32
cs_cin * cs_ci_nfree(cs_cin *N)
int * p
Definition: cs.h:68
cs_di * cs_di_spfree(cs_di *A)
cs_complex_t * x
Definition: evsl_cxsparse.c:34
double evsl_timer()
evsl timer for mac
Definition: mactime.c:14
int cs_di_lsolve(const cs_di *L, double *x)
cs_dis * cs_di_sfree(cs_dis *S)
cs_cis * cs_ci_sqr(int order, const cs_ci *A, int qr)