1 /*  sort.c  */
2 
3 #include "../SubMtx.h"
4 
5 #define MYDEBUG 0
6 
7 /*--------------------------------------------------------------------*/
8 /*
9    -----------------------------------------------------------
10    purpose -- sort the rows of the matrix into ascending order
11 
12    created -- 98mar02, cca
13    -----------------------------------------------------------
14 */
15 void
SubMtx_sortRowsUp(SubMtx * mtx)16 SubMtx_sortRowsUp (
17    SubMtx   *mtx
18 ) {
19 /*
20    ---------------
21    check the input
22    ---------------
23 */
24 if ( mtx == NULL ) {
25    fprintf(stderr, "\n fatal error in SubMtx_sortRowsUp(%p)"
26            "\n bad input\n", mtx) ;
27    exit(-1) ;
28 }
29 if ( ! (SUBMTX_IS_REAL(mtx) || SUBMTX_IS_COMPLEX(mtx)) ) {
30    fprintf(stderr, "\n fatal error in SubMtx_sortRowsUp(%p)"
31            "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
32            mtx, mtx->type) ;
33    exit(-1) ;
34 }
35 switch ( mtx->mode ) {
36 case SUBMTX_DENSE_ROWS :
37 case SUBMTX_DENSE_COLUMNS : {
38    A2       a2 ;
39    double   *entries ;
40    int      inc1, inc2, ncol, nrow ;
41    int      *rowind ;
42 
43    SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
44    A2_setDefaultFields(&a2) ;
45    A2_init(&a2, mtx->type, nrow, ncol, inc1, inc2, entries) ;
46    SubMtx_rowIndices(mtx, &nrow, &rowind) ;
47    A2_sortRowsUp(&a2, nrow, rowind) ;
48    } break ;
49 case SUBMTX_SPARSE_ROWS : {
50    double   *entries ;
51    int      ii, irow, jrow, kk, nent, nrow, offset, rowid, size ;
52    int      *indices, *ivtemp, *rowind, *sizes ;
53 
54    SubMtx_sparseRowsInfo(mtx, &nrow, &nent, &sizes, &indices, &entries);
55    SubMtx_rowIndices(mtx, &nrow, &rowind) ;
56 /*
57    ----------------------------------------------------------------
58    get a companion vector and fill with the row id's of the entries
59    ----------------------------------------------------------------
60 */
61    ivtemp = IVinit(nent, -1) ;
62    for ( irow = kk = 0 ; irow < nrow ; irow++ ) {
63       rowid = rowind[irow] ;
64       size  = sizes[irow] ;
65 #if MYDEBUG > 0
66       fprintf(stdout, "\n rowid %d, size %d, kk %d", rowid, size, kk) ;
67       fflush(stdout) ;
68 #endif
69       for ( ii = 0 ; ii < size ; ii++, kk++ ) {
70          ivtemp[kk] = rowid ;
71       }
72    }
73 #if MYDEBUG > 0
74    fprintf(stdout, "\n ivtemp[%d]", nent) ;
75    IVfprintf(stdout, nent, ivtemp) ;
76    fflush(stdout) ;
77 #endif
78 /*
79    -----------------------------------------------------------
80    zero the sizes vector, sort the (rowid,colid,entry) triples
81    into ascending order of rowids, and sort the row indices
82    -----------------------------------------------------------
83 */
84    IVzero(nrow, sizes) ;
85    if ( SUBMTX_IS_REAL(mtx) ) {
86       IV2DVqsortUp(nent, ivtemp, indices, entries) ;
87    } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
88       IV2ZVqsortUp(nent, ivtemp, indices, entries) ;
89    }
90 #if MYDEBUG > 0
91    fprintf(stdout, "\n after sort, ivtemp[%d]", nent) ;
92    IVfprintf(stdout, nent, ivtemp) ;
93    fflush(stdout) ;
94 #endif
95    IVqsortUp(nrow, rowind) ;
96 #if MYDEBUG > 0
97    fprintf(stdout, "\n after sort, rowind") ;
98    IVfprintf(stdout, nrow, rowind) ;
99    fflush(stdout) ;
100 #endif
101 /*
102    ----------------------------------------------------------------
103    sort each row in ascending order and fill the sizes[nrow] vector
104    ----------------------------------------------------------------
105 */
106    rowid = ivtemp[0] ;
107    jrow  = offset = 0 ;
108    kk    = size = 1 ;
109    while ( kk < nent ) {
110       if ( ivtemp[kk] == rowid ) {
111          size++ ;
112       } else {
113          while ( rowid != rowind[jrow] ) {
114 #if MYDEBUG > 0
115             fprintf(stdout, "\n rowid %d, rowind[%d] %d",
116                     rowid, jrow, rowind[jrow]) ;
117             fflush(stdout) ;
118 #endif
119             jrow++ ;
120          }
121          sizes[jrow] = size ;
122          if ( SUBMTX_IS_REAL(mtx) ) {
123             IVDVqsortUp(size, indices + offset, entries + offset) ;
124          } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
125             IVZVqsortUp(size, indices + offset, entries + 2*offset) ;
126          }
127          rowid = ivtemp[kk] ;
128          jrow++ ;
129          offset += size ;
130          size = 1 ;
131       }
132       kk++ ;
133    }
134    while ( rowid != rowind[jrow] ) {
135 #if MYDEBUG > 0
136       fprintf(stdout, "\n rowid %d, rowind[%d] %d",
137               rowid, jrow, rowind[jrow]) ;
138       fflush(stdout) ;
139 #endif
140       jrow++ ;
141    }
142    sizes[jrow] = size ;
143    if ( SUBMTX_IS_REAL(mtx) ) {
144       IVDVqsortUp(size, indices + offset, entries + offset) ;
145    } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
146       IVZVqsortUp(size, indices + offset, entries + 2*offset) ;
147    }
148    IVfree(ivtemp) ;
149    } break ;
150 default :
151    fprintf(stderr, "\n fatal error in SubMtx_sortRowsUp(%p)"
152            "\n bad type = %d", mtx, mtx->type) ;
153    exit(-1) ;
154 }
155 return ; }
156 
157 /*--------------------------------------------------------------------*/
158 /*
159    --------------------------------------------------------------
160    purpose -- sort the columns of the matrix into ascending order
161 
162    created -- 98mar02, cca
163    --------------------------------------------------------------
164 */
165 void
SubMtx_sortColumnsUp(SubMtx * mtx)166 SubMtx_sortColumnsUp (
167    SubMtx   *mtx
168 ) {
169 switch ( mtx->mode ) {
170 case SUBMTX_DENSE_ROWS :
171 case SUBMTX_DENSE_COLUMNS : {
172    A2       a2 ;
173    double   *entries ;
174    int      inc1, inc2, ncol, nrow ;
175    int      *colind ;
176 
177    A2_setDefaultFields(&a2) ;
178    SubMtx_denseInfo(mtx, &nrow, &ncol, &inc1, &inc2, &entries) ;
179    A2_init(&a2, mtx->type, nrow, ncol, inc1, inc2, entries) ;
180    SubMtx_columnIndices(mtx, &ncol, &colind) ;
181    A2_sortColumnsUp(&a2, ncol, colind) ;
182    } break ;
183 case SUBMTX_SPARSE_COLUMNS : {
184    double   *entries ;
185    int      colid, ii, jcol, kk, ncol, nent, offset, size ;
186    int      *colind, *indices, *ivtemp, *sizes ;
187 
188    SubMtx_sparseColumnsInfo(mtx,
189                           &ncol, &nent, &sizes, &indices, &entries) ;
190    SubMtx_columnIndices(mtx, &ncol, &colind) ;
191 /*
192    -------------------------------------------------------------------
193    get a companion vector and fill with the column id's of the entries
194    -------------------------------------------------------------------
195 */
196    ivtemp = IVinit(nent, -1) ;
197    for ( jcol = kk = 0 ; jcol < ncol ; jcol++ ) {
198       colid = colind[jcol] ;
199       size  = sizes[jcol] ;
200       for ( ii = 0 ; ii < size ; ii++, kk++ ) {
201          ivtemp[kk] = colid ;
202       }
203    }
204 #if MYDEBUG > 0
205    fprintf(stdout, "\n ivtemp[%d]", nent) ;
206    IVfprintf(stdout, nent, ivtemp) ;
207    fflush(stdout) ;
208 #endif
209 /*
210    -----------------------------------------------------------
211    zero the sizes vector, sort the (colid,rowid,entry) triples
212    into ascending order of colids, and sort the column indices
213    -----------------------------------------------------------
214 */
215    IVzero(ncol, sizes) ;
216    if ( SUBMTX_IS_REAL(mtx) ) {
217       IV2DVqsortUp(nent, ivtemp, indices, entries) ;
218    } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
219       IV2ZVqsortUp(nent, ivtemp, indices, entries) ;
220    }
221 #if MYDEBUG > 0
222    fprintf(stdout, "\n after sort, ivtemp[%d]", nent) ;
223    IVfprintf(stdout, nent, ivtemp) ;
224    fflush(stdout) ;
225 #endif
226    IVqsortUp(ncol, colind) ;
227 #if MYDEBUG > 0
228    fprintf(stdout, "\n after sort, colind") ;
229    IVfprintf(stdout, ncol, colind) ;
230    fflush(stdout) ;
231 #endif
232 /*
233    -------------------------------------------------------------------
234    sort each column in ascending order and fill the sizes[ncol] vector
235    -------------------------------------------------------------------
236 */
237    colid = ivtemp[0] ;
238    jcol  = offset = 0 ;
239    kk    = size = 1 ;
240    while ( kk < nent ) {
241       if ( ivtemp[kk] == colid ) {
242          size++ ;
243       } else {
244          while ( colid != colind[jcol] ) {
245 #if MYDEBUG > 0
246             fprintf(stdout, "\n colid %d, colind[%d] %d",
247                     colid, jcol, colind[jcol]) ;
248             fflush(stdout) ;
249 #endif
250             jcol++ ;
251          }
252          sizes[jcol] = size ;
253          if ( SUBMTX_IS_REAL(mtx) ) {
254             IVDVqsortUp(size, indices + offset, entries + offset) ;
255          } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
256             IVZVqsortUp(size, indices + offset, entries + 2*offset) ;
257          }
258          colid = ivtemp[kk] ;
259          jcol++ ;
260          offset += size ;
261          size = 1 ;
262       }
263       kk++ ;
264    }
265    while ( colid != colind[jcol] ) {
266 #if MYDEBUG > 0
267       fprintf(stdout, "\n colid %d, colind[%d] %d",
268               colid, jcol, colind[jcol]) ;
269       fflush(stdout) ;
270 #endif
271       jcol++ ;
272    }
273    sizes[jcol] = size ;
274    if ( SUBMTX_IS_REAL(mtx) ) {
275       IVDVqsortUp(size, indices + offset, entries + offset) ;
276    } else if ( SUBMTX_IS_COMPLEX(mtx) ) {
277       IVZVqsortUp(size, indices + offset, entries + 2*offset) ;
278    }
279    IVfree(ivtemp) ;
280    } break ;
281 default :
282    fprintf(stderr, "\n fatal error in SubMtx_sortColumnsUp(%p)"
283            "\n bad type = %d", mtx, mtx->type) ;
284    SubMtx_writeForHumanEye(mtx, stderr) ;
285    exit(-1) ;
286 }
287 return ; }
288 
289 /*--------------------------------------------------------------------*/
290