From c9cf2446c5eb18a29119c7dfe90c5b62f6b374c5 Mon Sep 17 00:00:00 2001 From: Justin Berger Date: Thu, 15 Mar 2018 23:14:29 -0600 Subject: Add sba poser --- redist/sba/sba_crsm.c | 338 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 338 insertions(+) create mode 100644 redist/sba/sba_crsm.c (limited to 'redist/sba/sba_crsm.c') 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 +#include + +#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(jcolidx[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(jcolidx[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(jcolidx[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; inr; ++i) + for(k=sm->rowptr[i]; krowptr[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; inr; ++i){ + low=sm->rowptr[i]; + high=sm->rowptr[i+1]-1; + if(jcolidx[low] || j>sm->colidx[high] || /* j not found */ + kcolidx[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(kcolidx[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