summaryrefslogtreecommitdiff
path: root/libmoped/libs/sba-1.6/sba_crsm.c
diff options
context:
space:
mode:
Diffstat (limited to 'libmoped/libs/sba-1.6/sba_crsm.c')
-rw-r--r--libmoped/libs/sba-1.6/sba_crsm.c346
1 files changed, 346 insertions, 0 deletions
diff --git a/libmoped/libs/sba-1.6/sba_crsm.c b/libmoped/libs/sba-1.6/sba_crsm.c
new file mode 100644
index 0000000..22179ea
--- /dev/null
+++ b/libmoped/libs/sba-1.6/sba_crsm.c
@@ -0,0 +1,346 @@
+/////////////////////////////////////////////////////////////////////////////////
+////
+//// 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