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