xref: /petsc/src/mat/impls/aij/seq/ij.c (revision 2d9e4a2ab1ebcf57c028f48d95136accae90829a)
1*2d9e4a2aSBarry Smith #ifndef lint
2*2d9e4a2aSBarry Smith static char vcid[] = "$Id: ij.c,v 1.1 1994/03/18 00:27:00 gropp Exp $";
3*2d9e4a2aSBarry Smith #endif
4*2d9e4a2aSBarry Smith 
5*2d9e4a2aSBarry Smith #include "aij.h"
6*2d9e4a2aSBarry Smith 
7*2d9e4a2aSBarry Smith /*
8*2d9e4a2aSBarry Smith   SpToSymmetricIJ - Convert a sparse AIJ matrix to IJ format
9*2d9e4a2aSBarry Smith            (ignore the "A" part)
10*2d9e4a2aSBarry Smith            Allocates the space needed. Uses only the lower triangular
11*2d9e4a2aSBarry Smith            part of the matrix.
12*2d9e4a2aSBarry Smith 
13*2d9e4a2aSBarry Smith     Description:
14*2d9e4a2aSBarry Smith     Take the data in the row-oriented sparse storage and build the
15*2d9e4a2aSBarry Smith     IJ data for the Matrix.  Return 0 on success, row + 1 on failure
16*2d9e4a2aSBarry Smith     at that row. Produces the ij for a symmetric matrix by only using
17*2d9e4a2aSBarry Smith     the lower triangular part of the matrix.
18*2d9e4a2aSBarry Smith 
19*2d9e4a2aSBarry Smith     Input Parameters:
20*2d9e4a2aSBarry Smith .   Matrix - matrix to convert
21*2d9e4a2aSBarry Smith 
22*2d9e4a2aSBarry Smith     Output Parameters:
23*2d9e4a2aSBarry Smith .   ia     - ia part of IJ representation (row information)
24*2d9e4a2aSBarry Smith .   ja     - ja part (column indices)
25*2d9e4a2aSBarry Smith 
26*2d9e4a2aSBarry Smith     Notes:
27*2d9e4a2aSBarry Smith     This routine is provided for ordering routines that require a
28*2d9e4a2aSBarry Smith     symmetric structure.  It is used in SpOrder (and derivatives) since
29*2d9e4a2aSBarry Smith     those routines call SparsePak routines that expect a symmetric
30*2d9e4a2aSBarry Smith     matrix.
31*2d9e4a2aSBarry Smith */
32*2d9e4a2aSBarry Smith int SpToSymmetricIJ( Matiaij *Matrix, int **iia, int **jja )
33*2d9e4a2aSBarry Smith {
34*2d9e4a2aSBarry Smith   int          *work,*ia,*ja,*j,i, nz, n, row, wr;
35*2d9e4a2aSBarry Smith   register int col;
36*2d9e4a2aSBarry Smith 
37*2d9e4a2aSBarry Smith   n  = Matrix->m;
38*2d9e4a2aSBarry Smith 
39*2d9e4a2aSBarry Smith   /* allocate space for row pointers */
40*2d9e4a2aSBarry Smith   *iia = ia = (int *) MALLOC( (n+1)*sizeof(int) ); CHKPTR(ia);
41*2d9e4a2aSBarry Smith   MEMSET(ia,0,(n+1)*sizeof(int));
42*2d9e4a2aSBarry Smith   work = (int *) MALLOC( (n+1)*sizeof(int) ); CHKPTR(work);
43*2d9e4a2aSBarry Smith 
44*2d9e4a2aSBarry Smith   /* determine the number of columns in each row */
45*2d9e4a2aSBarry Smith   ia[0] = 1;
46*2d9e4a2aSBarry Smith   for (row = 0; row < n; row++) {
47*2d9e4a2aSBarry Smith     nz = Matrix->i[row+1] - Matrix->i[row];
48*2d9e4a2aSBarry Smith     j  = Matrix->j + Matrix->i[row] - 1;
49*2d9e4a2aSBarry Smith     while (nz--) {
50*2d9e4a2aSBarry Smith        col = *j++;
51*2d9e4a2aSBarry Smith        if ( col > row ) {
52*2d9e4a2aSBarry Smith           break;
53*2d9e4a2aSBarry Smith        }
54*2d9e4a2aSBarry Smith        if (col != row) ia[row+1]++;
55*2d9e4a2aSBarry Smith        ia[col+1]++;
56*2d9e4a2aSBarry Smith     }
57*2d9e4a2aSBarry Smith  }
58*2d9e4a2aSBarry Smith 
59*2d9e4a2aSBarry Smith  /* shift ia[i] to point to next row */
60*2d9e4a2aSBarry Smith  for ( i=1; i<n+1; i++ ) {
61*2d9e4a2aSBarry Smith    row       = ia[i-1];
62*2d9e4a2aSBarry Smith    ia[i]     += row;
63*2d9e4a2aSBarry Smith    work[i-1] = row - 1;
64*2d9e4a2aSBarry Smith  }
65*2d9e4a2aSBarry Smith 
66*2d9e4a2aSBarry Smith  /* allocate space for column pointers */
67*2d9e4a2aSBarry Smith  nz = ia[n];
68*2d9e4a2aSBarry Smith  *jja = ja = (int *) MALLOC( nz*sizeof(int) ); CHKPTR(ja);
69*2d9e4a2aSBarry Smith 
70*2d9e4a2aSBarry Smith  /* loop over lower triangular part putting into ja */
71*2d9e4a2aSBarry Smith  for (row = 0; row < n; row++) {
72*2d9e4a2aSBarry Smith     nz = Matrix->i[row+1] - Matrix->i[row];
73*2d9e4a2aSBarry Smith     j  = Matrix->j + Matrix->i[row] - 1;
74*2d9e4a2aSBarry Smith     while (nz--) {
75*2d9e4a2aSBarry Smith        col = *j++;
76*2d9e4a2aSBarry Smith        if ( col > row ) {
77*2d9e4a2aSBarry Smith           break;
78*2d9e4a2aSBarry Smith        }
79*2d9e4a2aSBarry Smith        if (col != row) {wr = work[col]; work[col] = wr + 1;
80*2d9e4a2aSBarry Smith 			ja[wr] = row + 1; }
81*2d9e4a2aSBarry Smith        wr = work[row]; work[row] = wr + 1;
82*2d9e4a2aSBarry Smith        ja[wr] = col + 1;
83*2d9e4a2aSBarry Smith     }
84*2d9e4a2aSBarry Smith  }
85*2d9e4a2aSBarry Smith  FREE(work);
86*2d9e4a2aSBarry Smith   return 0;
87*2d9e4a2aSBarry Smith }
88*2d9e4a2aSBarry Smith 
89