EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
spmat.c
Go to the documentation of this file.
1 #include <string.h>
2 #include "def.h"
3 #include "blaslapack.h"
4 #include "struct.h"
5 #include "internal_proto.h"
6 
17 void csrcsc(int OUTINDEX, const int nrow, const int ncol, int job,
18  double *a, int *ja, int *ia,
19  double *ao, int *jao, int *iao) {
20  int i,k;
21  for (i=0; i<ncol+1; i++) {
22  iao[i] = 0;
23  }
24  // compute nnz of columns of A
25  for (i=0; i<nrow; i++) {
26  for (k=ia[i]; k<ia[i+1]; k++) {
27  iao[ja[k]+1] ++;
28  }
29  }
30  // compute pointers from lengths
31  for (i=0; i<ncol; i++) {
32  iao[i+1] += iao[i];
33  }
34  // now do the actual copying
35  for (i=0; i<nrow; i++) {
36  for (k=ia[i]; k<ia[i+1]; k++) {
37  int j = ja[k];
38  if (job) {
39  ao[iao[j]] = a[k];
40  }
41  jao[iao[j]++] = i + OUTINDEX;
42  }
43  }
44  /*---- reshift iao and leave */
45  for (i=ncol; i>0; i--) {
46  iao[i] = iao[i-1] + OUTINDEX;
47  }
48  iao[0] = OUTINDEX;
49 }
50 
56 void sortrow(csrMat *A) {
57  /*-------------------------------------------*/
58  int nrows = A->nrows;
59  int ncols = A->ncols;
60  int nnz = A->ia[nrows];
61  // work array
62  double *b;
63  int *jb, *ib;
64  Malloc(b, nnz, double);
65  Malloc(jb, nnz, int);
66  Malloc(ib, ncols+1, int);
67  // double transposition
68  csrcsc(0, nrows, ncols, 1, A->a, A->ja, A->ia, b, jb, ib);
69  csrcsc(0, ncols, nrows, 1, b, jb, ib, A->a, A->ja, A->ia);
70  // free
71  free(b);
72  free(jb);
73  free(ib);
74 }
75 
79 void csr_resize(int nrow, int ncol, int nnz, csrMat *csr) {
80  csr->owndata = 1;
81  csr->nrows = nrow;
82  csr->ncols = ncol;
83  Malloc(csr->ia, nrow+1, int);
84  Malloc(csr->ja, nnz, int);
85  Malloc(csr->a, nnz, double);
86 }
87 
91 void free_csr(csrMat *csr) {
92  /* if it does not own the data, do nothing */
93  if (!csr->owndata) {
94  return;
95  }
96  free(csr->ia);
97  free(csr->ja);
98  free(csr->a);
99 }
100 
106 void csr_copy(csrMat *A, csrMat *B, int allocB) {
107  int nrows = A->nrows;
108  int ncols = A->ncols;
109  int nnz = A->ia[nrows];
110  if (allocB) {
111  /* allocate memory */
112  csr_resize(nrows, ncols, nnz, B);
113  } else {
114  /* just set the sizes */
115  B->nrows = nrows;
116  B->ncols = ncols;
117  }
118  memcpy(B->ia, A->ia, (nrows+1)*sizeof(int));;
119  memcpy(B->ja, A->ja, nnz*sizeof(int));
120  memcpy(B->a, A->a, nnz*sizeof(double));
121 }
122 
126 void free_coo(cooMat *coo) {
127  free(coo->ir);
128  free(coo->jc);
129  free(coo->vv);
130 }
131 
135 int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr) {
136  const int nnz = coo->nnz;
137  //printf("@@@@ coo2csr, nnz %d\n", nnz);
138  /* allocate memory */
139  csr_resize(coo->nrows, coo->ncols, nnz, csr);
140  const int nrows = coo->nrows;
141  /* fill (ia, ja, a) */
142  int i;
143  for (i=0; i<nrows+1; i++) {
144  csr->ia[i] = 0;
145  }
146  for (i=0; i<nnz; i++) {
147  int row = coo->ir[i] - cooidx;
148  csr->ia[row+1] ++;
149  }
150  for (i=0; i<nrows; i++) {
151  csr->ia[i+1] += csr->ia[i];
152  }
153  for (i=0; i<nnz; i++) {
154  int row = coo->ir[i] - cooidx;
155  int col = coo->jc[i] - cooidx;
156  double val = coo->vv[i];
157  int k = csr->ia[row];
158  csr->a[k] = val;
159  csr->ja[k] = col;
160  csr->ia[row]++;
161  }
162  for (i=nrows; i>0; i--) {
163  csr->ia[i] = csr->ia[i-1];
164  }
165  csr->ia[0] = 0;
166  /* sort rows ? */
167  sortrow(csr);
168  return 0;
169 }
170 
171 
172 #if 0
173 
176 double dcsrinfnrm(csrMat *A){
177  // computes the inf-norm of A: max abs row sum
178  double ta = 0.0;
179  int nrows =A->nrows, i, j;
180  int *ia = A->ia;
181  double *aa = A->a;
182  /* for each row */
183  for (i=0; i<nrows; i++) {
184  /* abs row sum */
185  double t = 0.0;
186  for (j=ia[i]; j<ia[i+1]; j++) {
187  t += fabs(aa[j]);
188  }
189  /* take max */
190  ta = max(ta,t);
191  }
192  return (ta);
193 }
194 #endif
195 
199 void dcsrmv(char trans, int nrow, int ncol, double *a,
200  int *ia, int *ja, double *x, double *y) {
201  int len, jj=nrow;
202  if (trans == 'N') {
203  //#pragma omp parallel for schedule(guided)
204  double r;
205  /*for (i=0; i<nrow; i++) {
206  r = 0.0;
207  for (j=ia[i]; j<ia[i+1]; j++) {
208  r += a[j] * x[ja[j]];
209  }
210  y[i] = r;
211  }
212  */
213  while (jj--) {
214  len = *(ia+1) - *ia;
215  ia++;
216  r = 0.0;
217  while (len--)
218  r += x[*ja++]*(*a++);
219  *y++ = r;
220  }
221  } else {
222  double xi;
223  int jj, len;
224  jj = nrow;
225  //-------------------- this is from the matvec used in FILTLAN
226  // column oriented - gains up to 15% in time
227  double *w = y;
228  while (jj--)
229  *w++ = 0.0;
230  jj = nrow;
231  //-------------------- column loop
232  while (jj--){
233  len = *(ia+1) - *ia;
234  ia++;
235  xi = *x++;
236  while(len--)
237  y[*ja++] += xi*(*a++);
238  }
239  }
240 }
241 
247 void matvec_csr(double *x, double *y, void *data) {
248  csrMat *A = (csrMat *) data;
249 #ifdef USE_MKL
250  char cN = 'N';
251  /*
252  double alp = 1.0, bet = 0.0;
253  mkl_dcsrmv(&cN, &(A->nrows), &(A->ncols), &alp, "GXXCXX",
254  A->a, A->ja, A->ia, A->ia+1, x, &bet, y);
255  */
256  mkl_cspblas_dcsrgemv(&cN, &A->nrows, A->a, A->ia, A->ja, x, y);
257 #else
258  dcsrmv('N', A->nrows, A->ncols, A->a, A->ia, A->ja, x, y);
259 #endif
260 }
261 
262 
266 inline void matadd_insert(double t, csrMat *A, csrMat *C, int i, int *k,
267  int *j, int *map) {
268  /* if this entry already exists in C:
269  * checking if it is the first entry of this row
270  * and if column indices match
271  * NOTE that pointer j or k will be modified */
272  if (*k > C->ia[i] && C->ja[(*k)-1] == A->ja[*j]) {
273  if (map) {
274  /* j maps to k-1 in C */
275  map[(*j)] = (*k) - 1;
276  }
277  /* add to existing entry */
278  C->a[(*k)-1] += t * A->a[*j];
279  } else {
280  if (map) {
281  /* jA maps to k in C */
282  map[*j] = *k;
283  }
284  /* create new entry */
285  C->ja[*k] = A->ja[*j];
286  C->a[*k] = t * A->a[*j];
287  (*k)++;
288  }
289  (*j) ++;
290 }
291 
307 int matadd(double alp, double bet, csrMat *A, csrMat *B, csrMat *C,
308  int *mapA, int *mapB) {
309  int nnzA, nnzB, i, jA, jB, k;
310  /* check dimension */
311  if (A->nrows != B->nrows || A->ncols != B->ncols) {
312  return 1;
313  }
314  /* nnz of A and B */
315  nnzA = A->ia[A->nrows];
316  nnzB = B->ia[B->nrows];
317  /* alloc C [at most has nnzC = nnzA + nnzB] */
318  csr_resize(A->nrows, A->ncols, nnzA+nnzB, C);
319  /* nnz counter of C */
320  k = 0;
321  C->ia[0] = 0;
322  for (i=0; i<A->nrows; i++) {
323  /* open row i of A and B */
324  /* merging two sorted list */
325  for (jA=A->ia[i], jB=B->ia[i]; ; ) {
326  if (jA < A->ia[i+1] && jB < B->ia[i+1]) {
327  /* will insert the element with smaller col id */
328  if (A->ja[jA] <= B->ja[jB]) {
329  /* insert jA */
330  matadd_insert(alp, A, C, i, &k, &jA, mapA);
331  } else {
332  /* instert jB */
333  matadd_insert(bet, B, C, i, &k, &jB, mapB);
334  }
335  } else if (jA == A->ia[i+1]) {
336  for (; jB < B->ia[i+1]; ) {
337  /* instert jB */
338  matadd_insert(bet, B, C, i, &k, &jB, mapB);
339  }
340  break;
341  } else {
342  for (; jA < A->ia[i+1]; ) {
343  /* insert jA */
344  matadd_insert(alp, A, C, i, &k, &jA, mapA);
345  }
346  break;
347  }
348  }
349  C->ia[i+1] = k;
350  }
351  Realloc(C->ja, k, int);
352  Realloc(C->a, k, double);
353  return 0;
354 }
355 
357 int speye(int n, csrMat *A) {
358  int i;
359  csr_resize(n, n, n, A);
360  for (i=0; i<n; i++) {
361  A->ia[i] = A->ja[i] = i;
362  A->a[i] = 1.0;
363  }
364  A->ia[n] = n;
365  return 0;
366 }
367 
368 /*
369  * Diagonal scaling for A such that A = D^{-1}*A*D^{-1}, i.e.,
370  * A(i,j) = A(i,j) / (d(i)*d(j)), d = diag(D)
371  * @param[in,out] Coo Matrix to scale
372  * @param[in] d: The vector that contains d(i)
373  */
374 void diagScalCoo(cooMat *A, double *d) {
375  int i, row, col, nnz = A->nnz;
376  /* diagonal scaling for A */
377  for (i=0; i<nnz; i++) {
378  row = A->ir[i];
379  col = A->jc[i];
380  A->vv[i] /= d[row] * d[col];
381  }
382 }
383 
384 /*
385  * Diagonal scaling for A such that A = D^{-1}*A*D^{-1}, i.e.,
386  * A(i,j) = A(i,j) / (d(i)*d(j)), d = diag(D)
387  * @param[in,out] CSR Matrix to scale
388  * @param[in] d: The vector that contains d(i)
389  */
390 void diagScalCsr(csrMat *A, double *d) {
391  int i, j;
392  /* diagonal scaling for A */
393  for (i=0; i<A->nrows; i++) {
394  for (j=A->ia[i]; j<A->ia[i+1]; j++) {
395  A->a[j] /= d[i] * d[A->ja[j]];
396  }
397  }
398 }
399 
400 /*
401  * Extract the diagonal entries of csr matrix B
402  *
403  * @param[in] B Matrix to extract the diagonal
404  * @param[out] D preallocated vector of lengeth B.nrows
405  */
406 void extrDiagCsr(csrMat *B, double *d) {
407  int i, j, nrows = B->nrows;
408  for (i=0; i<nrows; i++) {
409  d[i] = 0.0;
410  for (j=B->ia[i]; j<B->ia[i+1]; j++) {
411  if (i == B->ja[j]) {
412  d[i] = B->a[j];
413  }
414  }
415  }
416 }
417 
418 /*
419  * Extract the upper triangular part of csr matrix A
420  *
421  * @param[in] A Matrix
422  * @param[out] U Matrix
423  */
424 void triuCsr(csrMat *A, csrMat *U) {
425  EVSL_Int i, j;
426  EVSL_Unsigned k = 0, nnzu = 0;
427 
428  /* count nnz in U */
429  for (i = 0; i < A->nrows; i++) {
430  for (j = A->ia[i]; j < A->ia[i+1]; j++) {
431  if (i <= A->ja[j]) {
432  nnzu++;
433  }
434  }
435  }
436  /* allocate memory for U */
437  csr_resize(A->nrows, A->ncols, nnzu, U);
438  U->ia[0] = 0;
439  /* copy entries from A to U */
440  for (i = 0; i < A->nrows; i++) {
441  for (j = A->ia[i]; j < A->ia[i+1]; j++) {
442  if (i <= A->ja[j]) {
443  U->ja[k] = A->ja[j];
444  U->a[k++] = A->a[j];
445  }
446  }
447  U->ia[i+1] = k;
448  }
449  CHKERR(k != nnzu);
450 }
451 
void dcsrmv(char trans, int nrow, int ncol, double *a, int *ia, int *ja, double *x, double *y)
csr matrix matvec or transpose matvec, (ia, ja, a) form
Definition: spmat.c:199
void triuCsr(csrMat *A, csrMat *U)
Definition: spmat.c:424
void csrcsc(int OUTINDEX, const int nrow, const int ncol, int job, double *a, int *ja, int *ia, double *ao, int *jao, int *iao)
convert csr to csc Assume input csr is 0-based index output csc 0/1 index specified by OUTINDEX * ...
Definition: spmat.c:17
void free_coo(cooMat *coo)
memory deallocation for coo matrix
Definition: spmat.c:126
void diagScalCoo(cooMat *A, double *d)
Definition: spmat.c:374
void extrDiagCsr(csrMat *B, double *d)
Definition: spmat.c:406
void csr_resize(int nrow, int ncol, int nnz, csrMat *csr)
memory allocation for csr matrix
Definition: spmat.c:79
defs in EVSL
int ncols
Definition: struct.h:32
#define EVSL_Unsigned
Definition: def.h:74
int nnz
Definition: struct.h:17
double * vv
Definition: struct.h:22
#define CHKERR(ierr)
Definition: def.h:19
int cooMat_to_csrMat(int cooidx, cooMat *coo, csrMat *csr)
convert coo to csr
Definition: spmat.c:135
void matadd_insert(double t, csrMat *A, csrMat *C, int i, int *k, int *j, int *map)
inline function used by matadd insert an element pointed by j of A (times t) to location k in C (row ...
Definition: spmat.c:266
int * jc
Definition: struct.h:17
void sortrow(csrMat *A)
Sort each row of a csr by increasing column order By double transposition.
Definition: spmat.c:56
int * ir
Definition: struct.h:17
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...
void diagScalCsr(csrMat *A, double *d)
Definition: spmat.c:390
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
structs used in evsl
sparse matrix format: the compressed sparse row (CSR) format, 0-based
Definition: struct.h:31
double * a
Definition: struct.h:37
int ncols
Definition: struct.h:17
#define max(a, b)
Definition: def.h:56
int * ia
Definition: struct.h:32
Defs for blaslapack routines.
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
int * ja
Definition: struct.h:32
int owndata
Definition: struct.h:32
int nrows
Definition: struct.h:17
void matvec_csr(double *x, double *y, void *data)
matvec for a CSR matrix, y = A*x. void *data points to csrMat, compatible form with EVSLMatvec (see s...
Definition: spmat.c:247
void csr_copy(csrMat *A, csrMat *B, int allocB)
copy a csr matrix A into B alloB: 0: will not allocate memory for B (have been alloced outside) 1: wi...
Definition: spmat.c:106
#define Realloc(base, nmem, type)
Definition: def.h:42
int nrows
Definition: struct.h:32
sparse matrix format: the coordinate (COO) format, 0-based
Definition: struct.h:16
#define EVSL_Int
Definition: def.h:73