1 /* ssrcsr.f -- translated by f2c (version 19970219).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 #include "f2c.h"
5
6 ssrcsr is part of the SPARSKIT library by Y. Saad at University of Minnesota
7 For the full SPARSKIT library visit:
8 www.cs.umn.edu/~saad/
9
10 */
11
ssrcsr(int * job,int * value2,int * nrow,double * a,int * ja,int * ia,int * nzmax,double * ao,int * jao,int * iao,int * indu,int * iwk,int * ierr)12 int ssrcsr(int *job, int *value2, int *nrow, double *a, int *ja, int *ia,
13 int *nzmax, double *ao, int *jao, int *iao, int *indu, int *iwk, int *ierr)
14 {
15 /* System generated locals */
16 int i__1, i__2, i__3;
17
18 /* Local variables */
19 static int ipos, i__, j, k, klast, kosav, ko, kfirst;
20 static double tmp;
21 static int nnz;
22
23 /* .. Scalar Arguments .. */
24 /* .. */
25 /* .. Array Arguments .. */
26 /* .. */
27 /* -----------------------------------------------------------------------
28 */
29 /* Symmetric Sparse Row to Compressed Sparse Row format */
30 /* -----------------------------------------------------------------------
31 */
32 /* This subroutine converts a given matrix in SSR format to regular */
33 /* CSR format by computing Ao = A + A' - diag(A), where A' is A */
34 /* transpose. */
35
36 /* Typically this routine is used to expand the SSR matrix of */
37 /* Harwell Boeing matrices, or to obtain a symmetrized graph of */
38 /* unsymmetric matrices. */
39
40 /* This routine is inplace, i.e., (Ao,jao,iao) may be same as */
41 /* (a,ja,ia). */
42
43 /* It is possible to input an arbitrary CSR matrix to this routine, */
44 /* since there is no syntactical difference between CSR and SSR */
45 /* format. It also removes duplicate entries and perform a partial */
46 /* ordering. The output matrix has an order of lower half, main */
47 /* diagonal and upper half after the partial ordering. */
48 /* -----------------------------------------------------------------------
49 */
50 /* on entry: */
51 /* --------- */
52
53 /* job = options */
54 /* 0 -- duplicate entries are not removed. If the input matrix is
55 */
56 /* SSR (not an arbitary CSR) matrix, no duplicate entry should
57 */
58 /* arise from this routine. */
59 /* 1 -- eliminate duplicate entries, zero entries. */
60 /* 2 -- eliminate duplicate entries and perform partial ordering.
61 */
62 /* 3 -- eliminate duplicate entries, sort the entries in the */
63 /* increasing order of clumn indices. */
64
65 /* value2= will the values of A be copied? */
66 /* 0 -- only expand the graph (a, ao are not touched) */
67 /* 1 -- expand the matrix with the values. */
68
69 /* nrow = column dimension of inout matrix */
70 /* a, */
71 /* ia, */
72 /* ja = matrix in compressed sparse row format. */
73
74 /* nzmax = size of arrays ao and jao. SSRCSR will abort if the storage */
75 /* provided in ao, jao is not sufficient to store A. See ierr. */
76
77 /* on return: */
78 /* ---------- */
79 /* ao, jao, iao */
80 /* = output matrix in compressed sparse row format. The resulting */
81 /* matrix is symmetric and is equal to A+A'-D. ao, jao, iao, */
82 /* can be the same as a, ja, ia in the calling sequence. */
83
84 /* indu = integer array of length nrow. INDU will contain pointers */
85 /* to the beginning of upper traigular part if job > 1. */
86 /* Otherwise it is also used as a work array (size nrow). */
87
88 /* iwk = integer work space (size nrow+1). */
89
90 /* ierr = integer. Serving as error message. If the length of the arrays
91 */
92 /* ao, jao exceeds nzmax, ierr returns the minimum value */
93 /* needed for nzmax. otherwise ierr=0 (normal return). */
94
95 /* -----------------------------------------------------------------------
96 */
97 /* .. Local Scalars .. */
98 /* .. */
99 /* .. Executable Statements .. */
100 /* Parameter adjustments */
101 --iwk;
102 --indu;
103 --iao;
104 --ia;
105 --a;
106 --ja;
107 --jao;
108 --ao;
109
110 /* Function Body */
111 *ierr = 0;
112 i__1 = *nrow;
113 for (i__ = 1; i__ <= i__1; ++i__) {
114 indu[i__] = 0;
115 iwk[i__] = 0;
116 /* L10: */
117 }
118 iwk[*nrow + 1] = 0;
119
120 /* .. compute number of elements in each row of (A'-D) */
121 /* put result in iwk(i+1) for row i. */
122
123 i__1 = *nrow;
124 for (i__ = 1; i__ <= i__1; ++i__) {
125 i__2 = ia[i__ + 1] - 1;
126 for (k = ia[i__]; k <= i__2; ++k) {
127 j = ja[k];
128 if (j != i__) {
129 ++iwk[j + 1];
130 }
131 /* L20: */
132 }
133 /* L30: */
134 }
135
136 /* .. find addresses of first elements of ouput matrix. result in iwk
137 */
138
139 iwk[1] = 1;
140 i__1 = *nrow;
141 for (i__ = 1; i__ <= i__1; ++i__) {
142 indu[i__] = iwk[i__] + ia[i__ + 1] - ia[i__];
143 iwk[i__ + 1] += indu[i__];
144 --indu[i__];
145 /* L40: */
146 }
147 /* .....Have we been given enough storage in ao, jao ? */
148 nnz = iwk[*nrow + 1] - 1;
149 if (nnz > *nzmax) {
150 *ierr = nnz;
151 return 0;
152 }
153
154 /* .. copy the existing matrix (backwards). */
155
156 kosav = iwk[*nrow + 1];
157 for (i__ = *nrow; i__ >= 1; --i__) {
158 klast = ia[i__ + 1] - 1;
159 kfirst = ia[i__];
160 iao[i__ + 1] = kosav;
161 kosav = iwk[i__];
162 ko = iwk[i__] - kfirst;
163 iwk[i__] = ko + klast + 1;
164 i__1 = kfirst;
165 for (k = klast; k >= i__1; --k) {
166 if (*value2 != 0) {
167 ao[k + ko] = a[k];
168 }
169 jao[k + ko] = ja[k];
170 /* L50: */
171 }
172 /* L60: */
173 }
174 iao[1] = 1;
175
176 /* now copy (A'-D). Go through the structure of ao, jao, iao */
177 /* that has already been copied. iwk(i) is the address */
178 /* of the next free location in row i for ao, jao. */
179
180 i__1 = *nrow;
181 for (i__ = 1; i__ <= i__1; ++i__) {
182 i__2 = indu[i__];
183 for (k = iao[i__]; k <= i__2; ++k) {
184 j = jao[k];
185 if (j != i__) {
186 ipos = iwk[j];
187 if (*value2 != 0) {
188 ao[ipos] = ao[k];
189 }
190 jao[ipos] = i__;
191 iwk[j] = ipos + 1;
192 }
193 /* L70: */
194 }
195 /* L80: */
196 }
197 if (*job <= 0) {
198 return 0;
199 }
200
201 /* .. eliminate duplicate entries -- */
202 /* array INDU is used as marker for existing indices, it is also the
203 */
204 /* location of the entry. */
205 /* IWK is used to stored the old IAO array. */
206 /* matrix is copied to squeeze out the space taken by the duplicated
207 */
208 /* entries. */
209
210 i__1 = *nrow;
211 for (i__ = 1; i__ <= i__1; ++i__) {
212 indu[i__] = 0;
213 iwk[i__] = iao[i__];
214 /* L90: */
215 }
216 iwk[*nrow + 1] = iao[*nrow + 1];
217 k = 1;
218 i__1 = *nrow;
219 for (i__ = 1; i__ <= i__1; ++i__) {
220 iao[i__] = k;
221 ipos = iwk[i__];
222 klast = iwk[i__ + 1];
223 L100:
224 if (ipos < klast) {
225 j = jao[ipos];
226 if (indu[j] == 0) {
227 /* .. new entry .. */
228 if (*value2 != 0) {
229 if (ao[ipos] != 0.) {
230 indu[j] = k;
231 jao[k] = jao[ipos];
232 ao[k] = ao[ipos];
233 ++k;
234 }
235 } else {
236 indu[j] = k;
237 jao[k] = jao[ipos];
238 ++k;
239 }
240 } else if (*value2 != 0) {
241 /* .. duplicate entry .. */
242 ao[indu[j]] += ao[ipos];
243 }
244 ++ipos;
245 goto L100;
246 }
247 /* .. remove marks before working on the next row .. */
248 i__2 = k - 1;
249 for (ipos = iao[i__]; ipos <= i__2; ++ipos) {
250 indu[jao[ipos]] = 0;
251 /* L110: */
252 }
253 /* L120: */
254 }
255 iao[*nrow + 1] = k;
256 if (*job <= 1) {
257 return 0;
258 }
259
260 /* .. partial ordering .. */
261 /* split the matrix into strict upper/lower triangular */
262 /* parts, INDU points to the the beginning of the strict upper part.
263 */
264
265 i__1 = *nrow;
266 for (i__ = 1; i__ <= i__1; ++i__) {
267 klast = iao[i__ + 1] - 1;
268 kfirst = iao[i__];
269 L130:
270 if (klast > kfirst) {
271 if (jao[klast] < i__ && jao[kfirst] >= i__) {
272 /* .. swap klast with kfirst .. */
273 j = jao[klast];
274 jao[klast] = jao[kfirst];
275 jao[kfirst] = j;
276 if (*value2 != 0) {
277 tmp = ao[klast];
278 ao[klast] = ao[kfirst];
279 ao[kfirst] = tmp;
280 }
281 }
282 if (jao[klast] >= i__) {
283 --klast;
284 }
285 if (jao[kfirst] < i__) {
286 ++kfirst;
287 }
288 goto L130;
289 }
290
291 if (jao[klast] < i__) {
292 indu[i__] = klast + 1;
293 } else {
294 indu[i__] = klast;
295 }
296 /* L140: */
297 }
298 if (*job <= 2) {
299 return 0;
300 }
301
302 /* .. order the entries according to column indices */
303 /* bubble-sort is used */
304
305 i__1 = *nrow;
306 for (i__ = 1; i__ <= i__1; ++i__) {
307 i__2 = indu[i__] - 1;
308 for (ipos = iao[i__]; ipos <= i__2; ++ipos) {
309 i__3 = ipos + 1;
310 for (j = indu[i__] - 1; j >= i__3; --j) {
311 k = j - 1;
312 if (jao[k] > jao[j]) {
313 ko = jao[k];
314 jao[k] = jao[j];
315 jao[j] = ko;
316 if (*value2 != 0) {
317 tmp = ao[k];
318 ao[k] = ao[j];
319 ao[j] = tmp;
320 }
321 }
322 /* L150: */
323 }
324 /* L160: */
325 }
326 i__2 = iao[i__ + 1] - 1;
327 for (ipos = indu[i__]; ipos <= i__2; ++ipos) {
328 i__3 = ipos + 1;
329 for (j = iao[i__ + 1] - 1; j >= i__3; --j) {
330 k = j - 1;
331 if (jao[k] > jao[j]) {
332 ko = jao[k];
333 jao[k] = jao[j];
334 jao[j] = ko;
335 if (*value2 != 0) {
336 tmp = ao[k];
337 ao[k] = ao[j];
338 ao[j] = tmp;
339 }
340 }
341 /* L170: */
342 }
343 /* L180: */
344 }
345 /* L190: */
346 }
347
348 return 0;
349 /* ---- end of ssrcsr ----------------------------------------------------
350 */
351 /* -----------------------------------------------------------------------
352 */
353 } /* ssrcsr_ */
354
355
356
357