aboutsummaryrefslogtreecommitdiff
path: root/redist/sba/sba_crsm.c
diff options
context:
space:
mode:
Diffstat (limited to 'redist/sba/sba_crsm.c')
-rw-r--r--redist/sba/sba_crsm.c338
1 files changed, 338 insertions, 0 deletions
diff --git a/redist/sba/sba_crsm.c b/redist/sba/sba_crsm.c
new file mode 100644
index 0000000..9ba11f1
--- /dev/null
+++ b/redist/sba/sba_crsm.c
@@ -0,0 +1,338 @@
+/////////////////////////////////////////////////////////////////////////////////
+////
+//// CRS sparse matrices manipulation routines
+//// Copyright (C) 2004-2008 Manolis Lourakis (lourakis at ics forth gr)
+//// Institute of Computer Science, Foundation for Research & Technology - Hellas
+//// Heraklion, Crete, Greece.
+////
+//// This program is free software; you can redistribute it and/or modify
+//// it under the terms of the GNU General Public License as published by
+//// the Free Software Foundation; either version 2 of the License, or
+//// (at your option) any later version.
+////
+//// This program is distributed in the hope that it will be useful,
+//// but WITHOUT ANY WARRANTY; without even the implied warranty of
+//// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+//// GNU General Public License for more details.
+////
+///////////////////////////////////////////////////////////////////////////////////
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "sba.h"
+
+static void sba_crsm_print(struct sba_crsm *sm, FILE *fp);
+static void sba_crsm_build(struct sba_crsm *sm, int *m, int nr, int nc);
+
+/* allocate a sparse CRS matrix */
+void sba_crsm_alloc(struct sba_crsm *sm, int nr, int nc, int nnz) {
+ int msz;
+
+ sm->nr = nr;
+ sm->nc = nc;
+ sm->nnz = nnz;
+ msz = 2 * nnz + nr + 1;
+ sm->val = (int *)malloc(msz * sizeof(int)); /* required memory is allocated in a single step */
+ if (!sm->val) {
+ fprintf(stderr, "SBA: memory allocation request failed in sba_crsm_alloc() [nr=%d, nc=%d, nnz=%d]\n", nr, nc,
+ nnz);
+ exit(1);
+ }
+ sm->colidx = sm->val + nnz;
+ sm->rowptr = sm->colidx + nnz;
+}
+
+/* free a sparse CRS matrix */
+void sba_crsm_free(struct sba_crsm *sm) {
+ sm->nr = sm->nc = sm->nnz = -1;
+ free(sm->val);
+ sm->val = sm->colidx = sm->rowptr = NULL;
+}
+
+static void sba_crsm_print(struct sba_crsm *sm, FILE *fp) {
+ register int i;
+
+ fprintf(fp, "matrix is %dx%d, %d non-zeros\nval: ", sm->nr, sm->nc, sm->nnz);
+ for (i = 0; i < sm->nnz; ++i)
+ fprintf(fp, "%d ", sm->val[i]);
+ fprintf(fp, "\ncolidx: ");
+ for (i = 0; i < sm->nnz; ++i)
+ fprintf(fp, "%d ", sm->colidx[i]);
+ fprintf(fp, "\nrowptr: ");
+ for (i = 0; i <= sm->nr; ++i)
+ fprintf(fp, "%d ", sm->rowptr[i]);
+ fprintf(fp, "\n");
+}
+
+/* build a sparse CRS matrix from a dense one. intended to serve as an example for sm creation */
+static void sba_crsm_build(struct sba_crsm *sm, int *m, int nr, int nc) {
+ int nnz;
+ register int i, j, k;
+
+ /* count nonzeros */
+ for (i = nnz = 0; i < nr; ++i)
+ for (j = 0; j < nc; ++j)
+ if (m[i * nc + j] != 0)
+ ++nnz;
+
+ sba_crsm_alloc(sm, nr, nc, nnz);
+
+ /* fill up the sm structure */
+ for (i = k = 0; i < nr; ++i) {
+ sm->rowptr[i] = k;
+ for (j = 0; j < nc; ++j)
+ if (m[i * nc + j] != 0) {
+ sm->val[k] = m[i * nc + j];
+ sm->colidx[k++] = j;
+ }
+ }
+ sm->rowptr[nr] = nnz;
+}
+
+/* returns the index of the (i, j) element. No bounds checking! */
+int sba_crsm_elmidx(struct sba_crsm *sm, int i, int j) {
+ register int low, high, mid, diff;
+
+ low = sm->rowptr[i];
+ high = sm->rowptr[i + 1] - 1;
+
+ /* binary search for finding the element at column j */
+ while (low <= high) {
+ /* following early termination test seems to actually slow down the search */
+ // if(j<sm->colidx[low] || j>sm->colidx[high]) return -1; /* not found */
+
+ /* mid=low+((high-low)>>1) ensures no index overflows */
+ mid = (low + high) >> 1; //(low+high)/2;
+ diff = j - sm->colidx[mid];
+ if (diff < 0)
+ high = mid - 1;
+ else if (diff > 0)
+ low = mid + 1;
+ else
+ return mid;
+ }
+
+ return -1; /* not found */
+}
+
+/* similarly to sba_crsm_elmidx() above, returns the index of the (i, j) element using the
+ * fact that the index of element (i, jp) was previously found to be jpidx. This can be
+ * slightly faster than sba_crsm_elmidx(). No bounds checking!
+ */
+int sba_crsm_elmidxp(struct sba_crsm *sm, int i, int j, int jp, int jpidx) {
+ register int low, high, mid, diff;
+
+ diff = j - jp;
+ if (diff > 0) {
+ low = jpidx + 1;
+ high = sm->rowptr[i + 1] - 1;
+ } else if (diff == 0)
+ return jpidx;
+ else { /* diff<0 */
+ low = sm->rowptr[i];
+ high = jpidx - 1;
+ }
+
+ /* binary search for finding the element at column j */
+ while (low <= high) {
+ /* following early termination test seems to actually slow down the search */
+ // if(j<sm->colidx[low] || j>sm->colidx[high]) return -1; /* not found */
+
+ /* mid=low+((high-low)>>1) ensures no index overflows */
+ mid = (low + high) >> 1; //(low+high)/2;
+ diff = j - sm->colidx[mid];
+ if (diff < 0)
+ high = mid - 1;
+ else if (diff > 0)
+ low = mid + 1;
+ else
+ return mid;
+ }
+
+ return -1; /* not found */
+}
+
+/* returns the number of nonzero elements in row i and
+ * fills up the vidxs and jidxs arrays with the val and column
+ * indexes of the elements found, respectively.
+ * vidxs and jidxs are assumed preallocated and of max. size sm->nc
+ */
+int sba_crsm_row_elmidxs(struct sba_crsm *sm, int i, int *vidxs, int *jidxs) {
+ register int j, k;
+
+ for (j = sm->rowptr[i], k = 0; j < sm->rowptr[i + 1]; ++j, ++k) {
+ vidxs[k] = j;
+ jidxs[k] = sm->colidx[j];
+ }
+
+ return k;
+}
+
+/* returns the number of nonzero elements in col j and
+ * fills up the vidxs and iidxs arrays with the val and row
+ * indexes of the elements found, respectively.
+ * vidxs and iidxs are assumed preallocated and of max. size sm->nr
+ */
+int sba_crsm_col_elmidxs(struct sba_crsm *sm, int j, int *vidxs, int *iidxs) {
+ register int *nextrowptr = sm->rowptr + 1;
+ register int i, l;
+ register int low, high, mid, diff;
+
+ for (i = l = 0; i < sm->nr; ++i) {
+ low = sm->rowptr[i];
+ high = nextrowptr[i] - 1;
+
+ /* binary search attempting to find an element at column j */
+ while (low <= high) {
+ // if(j<sm->colidx[low] || j>sm->colidx[high]) break; /* not found */
+
+ mid = (low + high) >> 1; //(low+high)/2;
+ diff = j - sm->colidx[mid];
+ if (diff < 0)
+ high = mid - 1;
+ else if (diff > 0)
+ low = mid + 1;
+ else { /* found */
+ vidxs[l] = mid;
+ iidxs[l++] = i;
+ break;
+ }
+ }
+ }
+
+ return l;
+}
+
+/* a more straighforward (but slower) implementation of the above function */
+/***
+int sba_crsm_col_elmidxs(struct sba_crsm *sm, int j, int *vidxs, int *iidxs)
+{
+register int i, k, l;
+
+ for(i=l=0; i<sm->nr; ++i)
+ for(k=sm->rowptr[i]; k<sm->rowptr[i+1]; ++k)
+ if(sm->colidx[k]==j){
+ vidxs[l]=k;
+ iidxs[l++]=i;
+ }
+
+ return l;
+}
+***/
+
+#if 0
+/* returns 1 if there exists a row i having columns j and k,
+ * i.e. a row i s.t. elements (i, j) and (i, k) are nonzero;
+ * 0 otherwise
+ */
+int sba_crsm_common_row(struct sba_crsm *sm, int j, int k)
+{
+register int i, low, high, mid, diff;
+
+ if(j==k) return 1;
+
+ for(i=0; i<sm->nr; ++i){
+ low=sm->rowptr[i];
+ high=sm->rowptr[i+1]-1;
+ if(j<sm->colidx[low] || j>sm->colidx[high] || /* j not found */
+ k<sm->colidx[low] || k>sm->colidx[high]) /* k not found */
+ continue;
+
+ /* binary search for finding the element at column j */
+ while(low<=high){
+ mid=(low+high)>>1; //(low+high)/2;
+ diff=j-sm->colidx[mid];
+ if(diff<0)
+ high=mid-1;
+ else if(diff>0)
+ low=mid+1;
+ else
+ goto jfound;
+ }
+
+ continue; /* j not found */
+
+jfound:
+ if(j>k){
+ low=sm->rowptr[i];
+ high=mid-1;
+ }
+ else{
+ low=mid+1;
+ high=sm->rowptr[i+1]-1;
+ }
+
+ if(k<sm->colidx[low] || k>sm->colidx[high]) continue; /* k not found */
+
+ /* binary search for finding the element at column k */
+ while(low<=high){
+ mid=(low+high)>>1; //(low+high)/2;
+ diff=k-sm->colidx[mid];
+ if(diff<0)
+ high=mid-1;
+ else if(diff>0)
+ low=mid+1;
+ else /* found */
+ return 1;
+ }
+ }
+
+ return 0;
+}
+#endif
+
+#if 0
+
+/* sample code using the above routines */
+
+main()
+{
+int mat[7][6]={
+ {10, 0, 0, 0, -2, 0},
+ {3, 9, 0, 0, 0, 3},
+ {0, 7, 8, 7, 0, 0},
+ {3, 0, 8, 7, 5, 0},
+ {0, 8, 0, 9, 9, 13},
+ {0, 4, 0, 0, 2, -1},
+ {3, 7, 0, 9, 2, 0}
+};
+
+struct sba_crsm sm;
+int i, j, k, l;
+int vidxs[7], /* max(6, 7) */
+ jidxs[6], iidxs[7];
+
+
+ sba_crsm_build(&sm, mat[0], 7, 6);
+ sba_crsm_print(&sm, stdout);
+
+ for(i=0; i<7; ++i){
+ for(j=0; j<6; ++j)
+ printf("%3d ", ((k=sba_crsm_elmidx(&sm, i, j))!=-1)? sm.val[k] : 0);
+ printf("\n");
+ }
+
+ for(i=0; i<7; ++i){
+ k=sba_crsm_row_elmidxs(&sm, i, vidxs, jidxs);
+ printf("row %d\n", i);
+ for(l=0; l<k; ++l){
+ j=jidxs[l];
+ printf("%d %d ", j, sm.val[vidxs[l]]);
+ }
+ printf("\n");
+ }
+
+ for(j=0; j<6; ++j){
+ k=sba_crsm_col_elmidxs(&sm, j, vidxs, iidxs);
+ printf("col %d\n", j);
+ for(l=0; l<k; ++l){
+ i=iidxs[l];
+ printf("%d %d ", i, sm.val[vidxs[l]]);
+ }
+ printf("\n");
+ }
+
+ sba_crsm_free(&sm);
+}
+#endif