1 /* sort.c */
2
3 #include "../A2.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 -----------------------------------------
8 permute the rows of the matrix
9 A(*,*) = A(index(*),*)
10 this method calls A2_sortRowsUp
11 but does not overwrite the index[] vector
12
13 created -- 98apr15, cca
14 -----------------------------------------
15 */
16 void
A2_permuteRows(A2 * mtx,int nrow,int index[])17 A2_permuteRows (
18 A2 *mtx,
19 int nrow,
20 int index[]
21 ) {
22 int *rowids ;
23 /*
24 ---------------
25 check the input
26 ---------------
27 */
28 if ( mtx == NULL || nrow < 0 || nrow > mtx->n1 || index == NULL ) {
29 fprintf(stderr, "\n fatal error in A2_permuteRows(%p,%d,%p)"
30 "\n bad input\n", mtx, nrow, index) ;
31 exit(-1) ;
32 }
33 rowids = IVinit(nrow, -1) ;
34 IVcopy(nrow, rowids, index) ;
35 A2_sortRowsUp(mtx, nrow, rowids) ;
36 IVfree(rowids) ;
37
38 return ; }
39
40 /*--------------------------------------------------------------------*/
41 /*
42 -----------------------------------------
43 permute the columns of the matrix
44 A(*,*) = A(*,index(*))
45 this method calls A2_sortColumnsUp
46 but does not overwrite the index[] vector
47
48 created -- 98apr15, cca
49 -----------------------------------------
50 */
51 void
A2_permuteColumns(A2 * mtx,int ncol,int index[])52 A2_permuteColumns (
53 A2 *mtx,
54 int ncol,
55 int index[]
56 ) {
57 int *colids ;
58 /*
59 ---------------
60 check the input
61 ---------------
62 */
63 if ( mtx == NULL || ncol < 0 || ncol > mtx->n2 || index == NULL ) {
64 fprintf(stderr, "\n fatal error in A2_permuteColumns(%p,%d,%p)"
65 "\n bad input\n", mtx, ncol, index) ;
66 exit(-1) ;
67 }
68 colids = IVinit(ncol, -1) ;
69 IVcopy(ncol, colids, index) ;
70 A2_sortColumnsUp(mtx, ncol, colids) ;
71 IVfree(colids) ;
72
73 return ; }
74
75 /*--------------------------------------------------------------------*/
76 /*
77 ----------------------------------------------
78 sort the rows of the matrix in ascending order
79 of the rowids[] vector. on return, rowids is
80 in asending order. return value is the number
81 of row swaps made.
82
83 created -- 98apr15, cca
84 ----------------------------------------------
85 */
86 int
A2_sortRowsUp(A2 * mtx,int nrow,int rowids[])87 A2_sortRowsUp (
88 A2 *mtx,
89 int nrow,
90 int rowids[]
91 ) {
92 int ii, minrow, minrowid, nswap, target ;
93 /*
94 ---------------
95 check the input
96 ---------------
97 */
98 if ( mtx == NULL || mtx->n1 < nrow || nrow < 0 || rowids == NULL ) {
99 fprintf(stderr, "\n fatal error in A2_sortRowsUp(%p,%d,%p)"
100 "\n bad input\n", mtx, nrow, rowids) ;
101 if ( mtx != NULL ) {
102 A2_writeStats(mtx, stderr) ;
103 }
104 exit(-1) ;
105 }
106 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
107 fprintf(stderr, "\n fatal error in A2_sortRowsUp(%p,%d,%p)"
108 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
109 mtx, nrow, rowids, mtx->type) ;
110 exit(-1) ;
111 }
112 nswap = 0 ;
113 if ( mtx->inc1 == 1 ) {
114 double *dvtmp ;
115 int jcol, ncol ;
116 int *ivtmp ;
117 /*
118 ---------------------------------------------------
119 matrix is stored by columns, so permute each column
120 ---------------------------------------------------
121 */
122 ivtmp = IVinit(nrow, -1) ;
123 if ( A2_IS_REAL(mtx) ) {
124 dvtmp = DVinit(nrow, 0.0) ;
125 } else if ( A2_IS_COMPLEX(mtx) ) {
126 dvtmp = DVinit(2*nrow, 0.0) ;
127 }
128 IVramp(nrow, ivtmp, 0, 1) ;
129 IV2qsortUp(nrow, rowids, ivtmp) ;
130 ncol = mtx->n2 ;
131 for ( jcol = 0 ; jcol < ncol ; jcol++ ) {
132 if ( A2_IS_REAL(mtx) ) {
133 DVcopy(nrow, dvtmp, A2_column(mtx, jcol)) ;
134 DVgather(nrow, A2_column(mtx, jcol), dvtmp, ivtmp) ;
135 } else if ( A2_IS_COMPLEX(mtx) ) {
136 ZVcopy(nrow, dvtmp, A2_column(mtx, jcol)) ;
137 ZVgather(nrow, A2_column(mtx, jcol), dvtmp, ivtmp) ;
138 }
139 }
140 IVfree(ivtmp) ;
141 DVfree(dvtmp) ;
142 } else {
143 /*
144 ----------------------------------------
145 use a simple insertion sort to swap rows
146 ----------------------------------------
147 */
148 for ( target = 0 ; target < nrow ; target++ ) {
149 minrow = target ;
150 minrowid = rowids[target] ;
151 for ( ii = target + 1 ; ii < nrow ; ii++ ) {
152 if ( minrowid > rowids[ii] ) {
153 minrow = ii ;
154 minrowid = rowids[ii] ;
155 }
156 }
157 if ( minrow != target ) {
158 rowids[minrow] = rowids[target] ;
159 rowids[target] = minrowid ;
160 A2_swapRows(mtx, target, minrow) ;
161 nswap++ ;
162 }
163 }
164 }
165
166 return(nswap) ; }
167
168 /*--------------------------------------------------------------------*/
169 /*
170 -------------------------------------------------
171 sort the columns of the matrix in ascending order
172 of the colids[] vector. on return, colids is
173 in asending order. return value is the number
174 of column swaps made.
175
176 created -- 98apr15, cca
177 -------------------------------------------------
178 */
179 int
A2_sortColumnsUp(A2 * mtx,int ncol,int colids[])180 A2_sortColumnsUp (
181 A2 *mtx,
182 int ncol,
183 int colids[]
184 ) {
185 int ii, mincol, mincolid, nswap, target ;
186 /*
187 ---------------
188 check the input
189 ---------------
190 */
191 if ( mtx == NULL || mtx->n2 < ncol || ncol < 0 || colids == NULL ) {
192 fprintf(stderr, "\n fatal error in A2_sortColumnsUp(%p,%d,%p)"
193 "\n bad input\n", mtx, ncol, colids) ;
194 if ( mtx != NULL ) {
195 A2_writeStats(mtx, stderr) ;
196 }
197 exit(-1) ;
198 }
199 if ( ! (A2_IS_REAL(mtx) || A2_IS_COMPLEX(mtx)) ) {
200 fprintf(stderr, "\n fatal error in A2_sortColumnsUp(%p,%d,%p)"
201 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
202 mtx, ncol, colids, mtx->type) ;
203 exit(-1) ;
204 }
205 nswap = 0 ;
206 if ( mtx->inc2 == 1 ) {
207 double *dvtmp ;
208 int irow, nrow ;
209 int *ivtmp ;
210 /*
211 ---------------------------------------------------
212 matrix is stored by rows, so permute each row
213 ---------------------------------------------------
214 */
215 ivtmp = IVinit(ncol, -1) ;
216 if ( A2_IS_REAL(mtx) ) {
217 dvtmp = DVinit(ncol, 0.0) ;
218 } else if ( A2_IS_COMPLEX(mtx) ) {
219 dvtmp = DVinit(2*ncol, 0.0) ;
220 }
221 IVramp(ncol, ivtmp, 0, 1) ;
222 IV2qsortUp(ncol, colids, ivtmp) ;
223 nrow = mtx->n1 ;
224 for ( irow = 0 ; irow < nrow ; irow++ ) {
225 if ( A2_IS_REAL(mtx) ) {
226 DVcopy(ncol, dvtmp, A2_row(mtx, irow)) ;
227 DVgather(ncol, A2_row(mtx, irow), dvtmp, ivtmp) ;
228 } else if ( A2_IS_COMPLEX(mtx) ) {
229 ZVcopy(ncol, dvtmp, A2_row(mtx, irow)) ;
230 ZVgather(ncol, A2_row(mtx, irow), dvtmp, ivtmp) ;
231 }
232 }
233 IVfree(ivtmp) ;
234 DVfree(dvtmp) ;
235 } else {
236 /*
237 ----------------------------------------
238 use a simple insertion sort to swap cols
239 ----------------------------------------
240 */
241 for ( target = 0 ; target < ncol ; target++ ) {
242 mincol = target ;
243 mincolid = colids[target] ;
244 for ( ii = target + 1 ; ii < ncol ; ii++ ) {
245 if ( mincolid > colids[ii] ) {
246 mincol = ii ;
247 mincolid = colids[ii] ;
248 }
249 }
250 if ( mincol != target ) {
251 colids[mincol] = colids[target] ;
252 colids[target] = mincolid ;
253 A2_swapColumns(mtx, target, mincol) ;
254 nswap++ ;
255 }
256 }
257 }
258 return(nswap) ; }
259
260 /*--------------------------------------------------------------------*/
261