Actual source code: unary.c
1: #define PETSCMAT_DLL
3: /* unary.f -- translated by f2c (version of 25 March 1992 12:58:56).
5: This code is protected by the GNU copyright. See the file
6: ilut.c in this directory for the full copyright. See below for the Author.
7: */
8: #include petsc.h
9: /* ----------------------------------------------------------------------- */
10: static PetscErrorCode SPARSEKIT2rperm(PetscInt *nrow,MatScalar *a,PetscInt *ja,PetscInt *ia,MatScalar *ao,PetscInt *jao,PetscInt *iao,PetscInt *perm,PetscInt *job)
11: {
12: /* System generated locals */
13: PetscInt i__1,i__2;
15: /* Local variables */
16: PetscInt i,j,k,ii,ko;
17: PetscInt values;
19: /* -----------------------------------------------------------------------
20: */
21: /* this subroutine permutes the rows of a matrix in CSR format. */
22: /* rperm computes B = P A where P is a permutation matrix. */
23: /* the permutation P is defined through the array perm: for each j, */
24: /* perm(j) represents the destination row number of row number j. */
25: /* Youcef Saad -- recoded Jan 28, 1991. */
26: /* -----------------------------------------------------------------------
27: */
28: /* on entry: */
29: /* ---------- */
30: /* n = dimension of the matrix */
31: /* a, ja, ia = input matrix in csr format */
32: /* perm = integer array of length nrow containing the permutation arrays
33: */
34: /* for the rows: perm(i) is the destination of row i in the */
35: /* permuted matrix. */
36: /* ---> a(i,j) in the original matrix becomes a(perm(i),j) */
37: /* in the output matrix. */
39: /* job = integer indicating the work to be done: */
40: /* job = 1 permute a, ja, ia into ao, jao, iao */
41: /* (including the copying of real values ao and */
42: /* the array iao). */
43: /* job .ne. 1 : ignore real values. */
44: /* (in which case arrays a and ao are not needed nor
45: */
46: /* used). */
48: /* ------------ */
49: /* on return: */
50: /* ------------ */
51: /* ao, jao, iao = input matrix in a, ja, ia format */
52: /* note : */
53: /* if (job.ne.1) then the arrays a and ao are not used. */
54: /* ----------------------------------------------------------------------c
55: */
56: /* Y. Saad, May 2, 1990 c
57: */
58: /* ----------------------------------------------------------------------c
59: */
60: /* Parameter adjustments */
61: --perm;
62: --iao;
63: --jao;
64: --ao;
65: --ia;
66: --ja;
67: --a;
69: /* Function Body */
70: values = *job == 1;
72: /* determine pointers for output matix. */
74: i__1 = *nrow;
75: for (j = 1; j <= i__1; ++j) {
76: i = perm[j];
77: iao[i + 1] = ia[j + 1] - ia[j];
78: /* L50: */
79: }
81: /* get pointers from lengths */
83: iao[1] = 1;
84: i__1 = *nrow;
85: for (j = 1; j <= i__1; ++j) {
86: iao[j + 1] += iao[j];
87: /* L51: */
88: }
90: /* copying */
92: i__1 = *nrow;
93: for (ii = 1; ii <= i__1; ++ii) {
95: /* old row = ii -- new row = iperm(ii) -- ko = new pointer */
97: ko = iao[perm[ii]];
98: i__2 = ia[ii + 1] - 1;
99: for (k = ia[ii]; k <= i__2; ++k) {
100: jao[ko] = ja[k];
101: if (values) {
102: ao[ko] = a[k];
103: }
104: ++ko;
105: /* L60: */
106: }
107: /* L100: */
108: }
110: return 0;
111: /* ---------end-of-rperm -------------------------------------------------
112: */
113: /* -----------------------------------------------------------------------
114: */
115: } /* rperm_ */
117: /* ----------------------------------------------------------------------- */
118: static PetscErrorCode SPARSEKIT2cperm(PetscInt *nrow,MatScalar * a,PetscInt * ja,PetscInt * ia,MatScalar * ao,PetscInt * jao,PetscInt * iao,PetscInt * perm,PetscInt * job)
119: {
120: /* System generated locals */
121: PetscInt i__1;
123: /* Local variables */
124: PetscInt i,k,nnz;
126: /* -----------------------------------------------------------------------
127: */
128: /* this subroutine permutes the columns of a matrix a, ja, ia. */
129: /* the result is written in the output matrix ao, jao, iao. */
130: /* cperm computes B = A P, where P is a permutation matrix */
131: /* that maps column j into column perm(j), i.e., on return */
132: /* a(i,j) becomes a(i,perm(j)) in new matrix */
133: /* Y. Saad, May 2, 1990 / modified Jan. 28, 1991. */
134: /* -----------------------------------------------------------------------
135: */
136: /* on entry: */
137: /* ---------- */
138: /* nrow = row dimension of the matrix */
140: /* a, ja, ia = input matrix in csr format. */
142: /* perm = integer array of length ncol (number of columns of A */
143: /* containing the permutation array the columns: */
144: /* a(i,j) in the original matrix becomes a(i,perm(j)) */
145: /* in the output matrix. */
147: /* job = integer indicating the work to be done: */
148: /* job = 1 permute a, ja, ia into ao, jao, iao */
149: /* (including the copying of real values ao and */
150: /* the array iao). */
151: /* job .ne. 1 : ignore real values ao and ignore iao. */
153: /* ------------ */
154: /* on return: */
155: /* ------------ */
156: /* ao, jao, iao = input matrix in a, ja, ia format (array ao not needed)
157: */
159: /* Notes: */
160: /* ------- */
161: /* 1. if job=1 then ao, iao are not used. */
162: /* 2. This routine is in place: ja, jao can be the same. */
163: /* 3. If the matrix is initially sorted (by increasing column number) */
164: /* then ao,jao,iao may not be on return. */
166: /* ----------------------------------------------------------------------c
167: */
168: /* local parameters: */
170: /* Parameter adjustments */
171: --perm;
172: --iao;
173: --jao;
174: --ao;
175: --ia;
176: --ja;
177: --a;
179: /* Function Body */
180: nnz = ia[*nrow + 1] - 1;
181: i__1 = nnz;
182: for (k = 1; k <= i__1; ++k) {
183: jao[k] = perm[ja[k]];
184: /* L100: */
185: }
187: /* done with ja array. return if no need to touch values. */
189: if (*job != 1) {
190: return 0;
191: }
193: /* else get new pointers -- and copy values too. */
195: i__1 = *nrow + 1;
196: for (i = 1; i <= i__1; ++i) {
197: iao[i] = ia[i];
198: /* L1: */
199: }
201: i__1 = nnz;
202: for (k = 1; k <= i__1; ++k) {
203: ao[k] = a[k];
204: /* L2: */
205: }
207: return 0;
208: /* ---------end-of-cperm--------------------------------------------------
209: */
210: /* -----------------------------------------------------------------------
211: */
212: } /* cperm_ */
214: /* ----------------------------------------------------------------------- */
215: PetscErrorCode SPARSEKIT2dperm(PetscInt *nrow,MatScalar *a,PetscInt *ja,PetscInt *ia,MatScalar *ao,PetscInt *jao,PetscInt *iao,PetscInt *perm,PetscInt *qperm,PetscInt *job)
216: {
217: PetscInt locjob;
219: /* -----------------------------------------------------------------------
220: */
221: /* This routine permutes the rows and columns of a matrix stored in CSR */
223: /* format. i.e., it computes P A Q, where P, Q are permutation matrices.
224: */
225: /* P maps row i into row perm(i) and Q maps column j into column qperm(j):
226: */
227: /* a(i,j) becomes a(perm(i),qperm(j)) in new matrix */
228: /* In the particular case where Q is the transpose of P (symmetric */
229: /* permutation of A) then qperm is not needed. */
230: /* note that qperm should be of length ncol (number of columns) but this
231: */
232: /* is not checked. */
233: /* -----------------------------------------------------------------------
234: */
235: /* Y. Saad, Sep. 21 1989 / recoded Jan. 28 1991. */
236: /* -----------------------------------------------------------------------
237: */
238: /* on entry: */
239: /* ---------- */
240: /* n = dimension of the matrix */
241: /* a, ja, */
242: /* ia = input matrix in a, ja, ia format */
243: /* perm = integer array of length n containing the permutation arrays */
244: /* for the rows: perm(i) is the destination of row i in the */
245: /* permuted matrix -- also the destination of column i in case */
246: /* permutation is symmetric (job .le. 2) */
248: /* qperm = same thing for the columns. This should be provided only */
249: /* if job=3 or job=4, i.e., only in the case of a nonsymmetric */
250: /* permutation of rows and columns. Otherwise qperm is a dummy */
252: /* job = integer indicating the work to be done: */
253: /* * job = 1,2 permutation is symmetric Ao :== P * A * transp(P) */
254: /* job = 1 permute a, ja, ia into ao, jao, iao */
255: /* job = 2 permute matrix ignoring real values. */
256: /* * job = 3,4 permutation is non-symmetric Ao :== P * A * Q */
257: /* job = 3 permute a, ja, ia into ao, jao, iao */
258: /* job = 4 permute matrix ignoring real values. */
260: /* on return: */
261: /* ----------- */
262: /* ao, jao, iao = input matrix in a, ja, ia format */
264: /* in case job .eq. 2 or job .eq. 4, a and ao are never referred to */
265: /* and can be dummy arguments. */
266: /* Notes: */
267: /* ------- */
268: /* 1) algorithm is in place */
269: /* 2) column indices may not be sorted on return even though they may be
270: */
271: /* on entry. */
272: /* ----------------------------------------------------------------------c
273: */
274: /* local variables */
276: /* locjob indicates whether or not real values must be copied. */
278: /* Parameter adjustments */
279: --qperm;
280: --perm;
281: --iao;
282: --jao;
283: --ao;
284: --ia;
285: --ja;
286: --a;
288: /* Function Body */
289: locjob = *job % 2;
291: /* permute rows first */
293: SPARSEKIT2rperm(nrow, &a[1], &ja[1], &ia[1], &ao[1], &jao[1], &iao[1], &perm[1], &locjob);
295: /* then permute columns */
297: locjob = 0;
299: if (*job <= 2) {
300: SPARSEKIT2cperm(nrow, &ao[1], &jao[1], &iao[1], &ao[1], &jao[1], &iao[1], &perm[1], &locjob);
301: } else {
302: SPARSEKIT2cperm(nrow, &ao[1], &jao[1], &iao[1], &ao[1], &jao[1], &iao[1], &qperm[1], &locjob);
303: }
305: return 0;
306: /* -------end-of-dperm----------------------------------------------------
307: */
308: /* -----------------------------------------------------------------------
309: */
310: } /* dperm_ */
312: /* ----------------------------------------------------------------------- */
313: PetscErrorCode SPARSEKIT2msrcsr(PetscInt *n,MatScalar * a,PetscInt * ja,MatScalar * ao,PetscInt * jao,PetscInt * iao,MatScalar * wk,PetscInt * iwk)
314: {
315: /* System generated locals */
316: PetscInt i__1, i__2;
318: /* Local variables */
319: PetscInt iptr;
320: PetscInt added;
321: PetscInt i, j, k, idiag, ii;
323: /* -----------------------------------------------------------------------
324: */
325: /* Modified - Sparse Row to Compressed Sparse Row */
327: /* -----------------------------------------------------------------------
328: */
329: /* converts a compressed matrix using a separated diagonal */
330: /* (modified sparse row format) in the Compressed Sparse Row */
331: /* format. */
332: /* does not check for zero elements in the diagonal. */
335: /* on entry : */
336: /* --------- */
337: /* n = row dimension of matrix */
338: /* a, ja = sparse matrix in msr sparse storage format */
339: /* see routine csrmsr for details on data structure */
341: /* on return : */
342: /* ----------- */
344: /* ao,jao,iao = output matrix in csr format. */
346: /* work arrays: */
347: /* ------------ */
348: /* wk = real work array of length n */
349: /* iwk = integer work array of length n+1 */
351: /* notes: */
352: /* The original version of this was NOT in place, but has */
353: /* been modified by adding the vector iwk to be in place. */
354: /* The original version had ja instead of iwk everywhere in */
355: /* loop 500. Modified Sun 29 May 1994 by R. Bramley (Indiana). */
357: /* -----------------------------------------------------------------------
358: */
359: /* Parameter adjustments */
360: --iwk;
361: --wk;
362: --iao;
363: --jao;
364: --ao;
365: --ja;
366: --a;
368: /* Function Body */
369: i__1 = *n;
370: for (i = 1; i <= i__1; ++i) {
371: wk[i] = a[i];
372: iwk[i] = ja[i];
373: /* L1: */
374: }
375: iwk[*n + 1] = ja[*n + 1];
376: iao[1] = 1;
377: iptr = 1;
378: /* --------- */
379: i__1 = *n;
380: for (ii = 1; ii <= i__1; ++ii) {
381: added = 0;
382: idiag = iptr + (iwk[ii + 1] - iwk[ii]);
383: i__2 = iwk[ii + 1] - 1;
384: for (k = iwk[ii]; k <= i__2; ++k) {
385: j = ja[k];
386: if (j < ii) {
387: ao[iptr] = a[k];
388: jao[iptr] = j;
389: ++iptr;
390: } else if (added) {
391: ao[iptr] = a[k];
392: jao[iptr] = j;
393: ++iptr;
394: } else {
395: /* add diag element - only reserve a position for it. */
396: idiag = iptr;
397: ++iptr;
398: added = 1;
399: /* then other element */
400: ao[iptr] = a[k];
401: jao[iptr] = j;
402: ++iptr;
403: }
404: /* L100: */
405: }
406: ao[idiag] = wk[ii];
407: jao[idiag] = ii;
408: if (! added) {
409: ++iptr;
410: }
411: iao[ii + 1] = iptr;
412: /* L500: */
413: }
414: return 0;
415: } /* msrcsr_ */