Warning fixes, one actual bug found in sequencer sound wave drawing. Also
[blender-staging.git] / intern / opennl / superlu / get_perm_c.c
1 /*
2  * -- SuperLU routine (version 2.0) --
3  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
4  * and Lawrence Berkeley National Lab.
5  * November 15, 1997
6  *
7  */
8
9 #include "ssp_defs.h"
10 #include "colamd.h"
11
12 extern int  genmmd_(int *, int *, int *, int *, int *, int *, int *, 
13                     int *, int *, int *, int *, int *);
14
15 static void
16 get_colamd(
17            const int m,  /* number of rows in matrix A. */
18            const int n,  /* number of columns in matrix A. */
19            const int nnz,/* number of nonzeros in matrix A. */
20            int *colptr,  /* column pointer of size n+1 for matrix A. */
21            int *rowind,  /* row indices of size nz for matrix A. */
22            int *perm_c   /* out - the column permutation vector. */
23            )
24 {
25     int Alen, *A, i, info, *p;
26     double *knobs;
27
28     Alen = colamd_recommended(nnz, m, n);
29
30     if ( !(knobs = (double *) SUPERLU_MALLOC(COLAMD_KNOBS * sizeof(double))) )
31         ABORT("Malloc fails for knobs");
32     colamd_set_defaults(knobs);
33
34     if (!(A = (int *) SUPERLU_MALLOC(Alen * sizeof(int))) )
35         ABORT("Malloc fails for A[]");
36     if (!(p = (int *) SUPERLU_MALLOC((n+1) * sizeof(int))) )
37         ABORT("Malloc fails for p[]");
38     for (i = 0; i <= n; ++i) p[i] = colptr[i];
39     for (i = 0; i < nnz; ++i) A[i] = rowind[i];
40     info = colamd(m, n, Alen, A, p, knobs);
41     if ( info == FALSE ) ABORT("COLAMD failed");
42
43     for (i = 0; i < n; ++i) perm_c[p[i]] = i;
44
45     SUPERLU_FREE(knobs);
46     SUPERLU_FREE(A);
47     SUPERLU_FREE(p);
48 }
49
50 static void
51 getata(
52        const int m,      /* number of rows in matrix A. */
53        const int n,      /* number of columns in matrix A. */
54        const int nz,     /* number of nonzeros in matrix A */
55        int *colptr,      /* column pointer of size n+1 for matrix A. */
56        int *rowind,      /* row indices of size nz for matrix A. */
57        int *atanz,       /* out - on exit, returns the actual number of
58                             nonzeros in matrix A'*A. */
59        int **ata_colptr, /* out - size n+1 */
60        int **ata_rowind  /* out - size *atanz */
61        )
62 /*
63  * Purpose
64  * =======
65  *
66  * Form the structure of A'*A. A is an m-by-n matrix in column oriented
67  * format represented by (colptr, rowind). The output A'*A is in column
68  * oriented format (symmetrically, also row oriented), represented by
69  * (ata_colptr, ata_rowind).
70  *
71  * This routine is modified from GETATA routine by Tim Davis.
72  * The complexity of this algorithm is: SUM_{i=1,m} r(i)^2,
73  * i.e., the sum of the square of the row counts.
74  *
75  * Questions
76  * =========
77  *     o  Do I need to withhold the *dense* rows?
78  *     o  How do I know the number of nonzeros in A'*A?
79  * 
80  */
81 {
82     register int i, j, k, col, num_nz, ti, trow;
83     int *marker, *b_colptr, *b_rowind;
84     int *t_colptr, *t_rowind; /* a column oriented form of T = A' */
85
86     if ( !(marker = (int*) SUPERLU_MALLOC((SUPERLU_MAX(m,n)+1)*sizeof(int))) )
87         ABORT("SUPERLU_MALLOC fails for marker[]");
88     if ( !(t_colptr = (int*) SUPERLU_MALLOC((m+1) * sizeof(int))) )
89         ABORT("SUPERLU_MALLOC t_colptr[]");
90     if ( !(t_rowind = (int*) SUPERLU_MALLOC(nz * sizeof(int))) )
91         ABORT("SUPERLU_MALLOC fails for t_rowind[]");
92
93     
94     /* Get counts of each column of T, and set up column pointers */
95     for (i = 0; i < m; ++i) marker[i] = 0;
96     for (j = 0; j < n; ++j) {
97         for (i = colptr[j]; i < colptr[j+1]; ++i)
98             ++marker[rowind[i]];
99     }
100     t_colptr[0] = 0;
101     for (i = 0; i < m; ++i) {
102         t_colptr[i+1] = t_colptr[i] + marker[i];
103         marker[i] = t_colptr[i];
104     }
105
106     /* Transpose the matrix from A to T */
107     for (j = 0; j < n; ++j)
108         for (i = colptr[j]; i < colptr[j+1]; ++i) {
109             col = rowind[i];
110             t_rowind[marker[col]] = j;
111             ++marker[col];
112         }
113
114     
115     /* ----------------------------------------------------------------
116        compute B = T * A, where column j of B is:
117
118        Struct (B_*j) =    UNION   ( Struct (T_*k) )
119                         A_kj != 0
120
121        do not include the diagonal entry
122    
123        ( Partition A as: A = (A_*1, ..., A_*n)
124          Then B = T * A = (T * A_*1, ..., T * A_*n), where
125          T * A_*j = (T_*1, ..., T_*m) * A_*j.  )
126        ---------------------------------------------------------------- */
127
128     /* Zero the diagonal flag */
129     for (i = 0; i < n; ++i) marker[i] = -1;
130
131     /* First pass determines number of nonzeros in B */
132     num_nz = 0;
133     for (j = 0; j < n; ++j) {
134         /* Flag the diagonal so it's not included in the B matrix */
135         marker[j] = j;
136
137         for (i = colptr[j]; i < colptr[j+1]; ++i) {
138             /* A_kj is nonzero, add pattern of column T_*k to B_*j */
139             k = rowind[i];
140             for (ti = t_colptr[k]; ti < t_colptr[k+1]; ++ti) {
141                 trow = t_rowind[ti];
142                 if ( marker[trow] != j ) {
143                     marker[trow] = j;
144                     num_nz++;
145                 }
146             }
147         }
148     }
149     *atanz = num_nz;
150     
151     /* Allocate storage for A'*A */
152     if ( !(*ata_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
153         ABORT("SUPERLU_MALLOC fails for ata_colptr[]");
154     if ( *atanz ) {
155         if ( !(*ata_rowind = (int*) SUPERLU_MALLOC( *atanz * sizeof(int)) ) )
156             ABORT("SUPERLU_MALLOC fails for ata_rowind[]");
157     }
158     b_colptr = *ata_colptr; /* aliasing */
159     b_rowind = *ata_rowind;
160     
161     /* Zero the diagonal flag */
162     for (i = 0; i < n; ++i) marker[i] = -1;
163     
164     /* Compute each column of B, one at a time */
165     num_nz = 0;
166     for (j = 0; j < n; ++j) {
167         b_colptr[j] = num_nz;
168         
169         /* Flag the diagonal so it's not included in the B matrix */
170         marker[j] = j;
171
172         for (i = colptr[j]; i < colptr[j+1]; ++i) {
173             /* A_kj is nonzero, add pattern of column T_*k to B_*j */
174             k = rowind[i];
175             for (ti = t_colptr[k]; ti < t_colptr[k+1]; ++ti) {
176                 trow = t_rowind[ti];
177                 if ( marker[trow] != j ) {
178                     marker[trow] = j;
179                     b_rowind[num_nz++] = trow;
180                 }
181             }
182         }
183     }
184     b_colptr[n] = num_nz;
185        
186     SUPERLU_FREE(marker);
187     SUPERLU_FREE(t_colptr);
188     SUPERLU_FREE(t_rowind);
189 }
190
191
192 static void
193 at_plus_a(
194           const int n,      /* number of columns in matrix A. */
195           const int nz,     /* number of nonzeros in matrix A */
196           int *colptr,      /* column pointer of size n+1 for matrix A. */
197           int *rowind,      /* row indices of size nz for matrix A. */
198           int *bnz,         /* out - on exit, returns the actual number of
199                                nonzeros in matrix A'*A. */
200           int **b_colptr,   /* out - size n+1 */
201           int **b_rowind    /* out - size *bnz */
202           )
203 {
204 /*
205  * Purpose
206  * =======
207  *
208  * Form the structure of A'+A. A is an n-by-n matrix in column oriented
209  * format represented by (colptr, rowind). The output A'+A is in column
210  * oriented format (symmetrically, also row oriented), represented by
211  * (b_colptr, b_rowind).
212  *
213  */
214     register int i, j, k, col, num_nz;
215     int *t_colptr, *t_rowind; /* a column oriented form of T = A' */
216     int *marker;
217
218     if ( !(marker = (int*) SUPERLU_MALLOC( n * sizeof(int)) ) )
219         ABORT("SUPERLU_MALLOC fails for marker[]");
220     if ( !(t_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
221         ABORT("SUPERLU_MALLOC fails for t_colptr[]");
222     if ( !(t_rowind = (int*) SUPERLU_MALLOC( nz * sizeof(int)) ) )
223         ABORT("SUPERLU_MALLOC fails t_rowind[]");
224
225     
226     /* Get counts of each column of T, and set up column pointers */
227     for (i = 0; i < n; ++i) marker[i] = 0;
228     for (j = 0; j < n; ++j) {
229         for (i = colptr[j]; i < colptr[j+1]; ++i)
230             ++marker[rowind[i]];
231     }
232     t_colptr[0] = 0;
233     for (i = 0; i < n; ++i) {
234         t_colptr[i+1] = t_colptr[i] + marker[i];
235         marker[i] = t_colptr[i];
236     }
237
238     /* Transpose the matrix from A to T */
239     for (j = 0; j < n; ++j)
240         for (i = colptr[j]; i < colptr[j+1]; ++i) {
241             col = rowind[i];
242             t_rowind[marker[col]] = j;
243             ++marker[col];
244         }
245
246
247     /* ----------------------------------------------------------------
248        compute B = A + T, where column j of B is:
249
250        Struct (B_*j) = Struct (A_*k) UNION Struct (T_*k)
251
252        do not include the diagonal entry
253        ---------------------------------------------------------------- */
254
255     /* Zero the diagonal flag */
256     for (i = 0; i < n; ++i) marker[i] = -1;
257
258     /* First pass determines number of nonzeros in B */
259     num_nz = 0;
260     for (j = 0; j < n; ++j) {
261         /* Flag the diagonal so it's not included in the B matrix */
262         marker[j] = j;
263
264         /* Add pattern of column A_*k to B_*j */
265         for (i = colptr[j]; i < colptr[j+1]; ++i) {
266             k = rowind[i];
267             if ( marker[k] != j ) {
268                 marker[k] = j;
269                 ++num_nz;
270             }
271         }
272
273         /* Add pattern of column T_*k to B_*j */
274         for (i = t_colptr[j]; i < t_colptr[j+1]; ++i) {
275             k = t_rowind[i];
276             if ( marker[k] != j ) {
277                 marker[k] = j;
278                 ++num_nz;
279             }
280         }
281     }
282     *bnz = num_nz;
283     
284     /* Allocate storage for A+A' */
285     if ( !(*b_colptr = (int*) SUPERLU_MALLOC( (n+1) * sizeof(int)) ) )
286         ABORT("SUPERLU_MALLOC fails for b_colptr[]");
287     if ( *bnz) {
288       if ( !(*b_rowind = (int*) SUPERLU_MALLOC( *bnz * sizeof(int)) ) )
289         ABORT("SUPERLU_MALLOC fails for b_rowind[]");
290     }
291     
292     /* Zero the diagonal flag */
293     for (i = 0; i < n; ++i) marker[i] = -1;
294     
295     /* Compute each column of B, one at a time */
296     num_nz = 0;
297     for (j = 0; j < n; ++j) {
298         (*b_colptr)[j] = num_nz;
299         
300         /* Flag the diagonal so it's not included in the B matrix */
301         marker[j] = j;
302
303         /* Add pattern of column A_*k to B_*j */
304         for (i = colptr[j]; i < colptr[j+1]; ++i) {
305             k = rowind[i];
306             if ( marker[k] != j ) {
307                 marker[k] = j;
308                 (*b_rowind)[num_nz++] = k;
309             }
310         }
311
312         /* Add pattern of column T_*k to B_*j */
313         for (i = t_colptr[j]; i < t_colptr[j+1]; ++i) {
314             k = t_rowind[i];
315             if ( marker[k] != j ) {
316                 marker[k] = j;
317                 (*b_rowind)[num_nz++] = k;
318             }
319         }
320     }
321     (*b_colptr)[n] = num_nz;
322        
323     SUPERLU_FREE(marker);
324     SUPERLU_FREE(t_colptr);
325     SUPERLU_FREE(t_rowind);
326 }
327
328 void
329 get_perm_c(int ispec, SuperMatrix *A, int *perm_c)
330 /*
331  * Purpose
332  * =======
333  *
334  * GET_PERM_C obtains a permutation matrix Pc, by applying the multiple
335  * minimum degree ordering code by Joseph Liu to matrix A'*A or A+A'.
336  * or using approximate minimum degree column ordering by Davis et. al.
337  * The LU factorization of A*Pc tends to have less fill than the LU 
338  * factorization of A.
339  *
340  * Arguments
341  * =========
342  *
343  * ispec   (input) int
344  *         Specifies the type of column ordering to reduce fill:
345  *         = 1: minimum degree on the structure of A^T * A
346  *         = 2: minimum degree on the structure of A^T + A
347  *         = 3: approximate minimum degree for unsymmetric matrices
348  *         If ispec == 0, the natural ordering (i.e., Pc = I) is returned.
349  * 
350  * A       (input) SuperMatrix*
351  *         Matrix A in A*X=B, of dimension (A->nrow, A->ncol). The number
352  *         of the linear equations is A->nrow. Currently, the type of A 
353  *         can be: Stype = NC; Dtype = _D; Mtype = GE. In the future,
354  *         more general A can be handled.
355  *
356  * perm_c  (output) int*
357  *         Column permutation vector of size A->ncol, which defines the 
358  *         permutation matrix Pc; perm_c[i] = j means column i of A is 
359  *         in position j in A*Pc.
360  *
361  */
362 {
363     NCformat *Astore = A->Store;
364     int m, n, bnz, *b_colptr, i;
365     int delta, maxint, nofsub, *invp;
366     int *b_rowind, *dhead, *qsize, *llist, *marker;
367     double t, SuperLU_timer_();
368     
369     /* make gcc happy */
370     b_rowind=NULL;
371     b_colptr=NULL;
372     
373     m = A->nrow;
374     n = A->ncol;
375
376     t = SuperLU_timer_();
377     switch ( ispec ) {
378         case 0: /* Natural ordering */
379               for (i = 0; i < n; ++i) perm_c[i] = i;
380 #if ( PRNTlevel>=1 )
381               printf("Use natural column ordering.\n");
382 #endif
383               return;
384         case 1: /* Minimum degree ordering on A'*A */
385               getata(m, n, Astore->nnz, Astore->colptr, Astore->rowind,
386                      &bnz, &b_colptr, &b_rowind);
387 #if ( PRNTlevel>=1 )
388               printf("Use minimum degree ordering on A'*A.\n");
389 #endif
390               t = SuperLU_timer_() - t;
391               /*printf("Form A'*A time = %8.3f\n", t);*/
392               break;
393         case 2: /* Minimum degree ordering on A'+A */
394               if ( m != n ) ABORT("Matrix is not square");
395               at_plus_a(n, Astore->nnz, Astore->colptr, Astore->rowind,
396                         &bnz, &b_colptr, &b_rowind);
397 #if ( PRNTlevel>=1 )
398               printf("Use minimum degree ordering on A'+A.\n");
399 #endif
400               t = SuperLU_timer_() - t;
401               /*printf("Form A'+A time = %8.3f\n", t);*/
402               break;
403         case 3: /* Approximate minimum degree column ordering. */
404               get_colamd(m, n, Astore->nnz, Astore->colptr, Astore->rowind,
405                          perm_c);
406 #if ( PRNTlevel>=1 )
407               printf(".. Use approximate minimum degree column ordering.\n");
408 #endif
409               return; 
410         default:
411               ABORT("Invalid ISPEC");
412                   return;
413     }
414
415     if ( bnz != 0 ) {
416         t = SuperLU_timer_();
417
418         /* Initialize and allocate storage for GENMMD. */
419         delta = 1; /* DELTA is a parameter to allow the choice of nodes
420                       whose degree <= min-degree + DELTA. */
421         maxint = 2147483647; /* 2**31 - 1 */
422         invp = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
423         if ( !invp ) ABORT("SUPERLU_MALLOC fails for invp.");
424         dhead = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
425         if ( !dhead ) ABORT("SUPERLU_MALLOC fails for dhead.");
426         qsize = (int *) SUPERLU_MALLOC((n+delta)*sizeof(int));
427         if ( !qsize ) ABORT("SUPERLU_MALLOC fails for qsize.");
428         llist = (int *) SUPERLU_MALLOC(n*sizeof(int));
429         if ( !llist ) ABORT("SUPERLU_MALLOC fails for llist.");
430         marker = (int *) SUPERLU_MALLOC(n*sizeof(int));
431         if ( !marker ) ABORT("SUPERLU_MALLOC fails for marker.");
432
433         /* Transform adjacency list into 1-based indexing required by GENMMD.*/
434         for (i = 0; i <= n; ++i) ++b_colptr[i];
435         for (i = 0; i < bnz; ++i) ++b_rowind[i];
436         
437         genmmd_(&n, b_colptr, b_rowind, perm_c, invp, &delta, dhead, 
438                 qsize, llist, marker, &maxint, &nofsub);
439
440         /* Transform perm_c into 0-based indexing. */
441         for (i = 0; i < n; ++i) --perm_c[i];
442
443         SUPERLU_FREE(b_colptr);
444         SUPERLU_FREE(b_rowind);
445         SUPERLU_FREE(invp);
446         SUPERLU_FREE(dhead);
447         SUPERLU_FREE(qsize);
448         SUPERLU_FREE(llist);
449         SUPERLU_FREE(marker);
450
451         t = SuperLU_timer_() - t;
452         /*  printf("call GENMMD time = %8.3f\n", t);*/
453
454     } else { /* Empty adjacency structure */
455         for (i = 0; i < n; ++i) perm_c[i] = i;
456     }
457
458 }