32455dd517148ab1ae8768fb20f813b6ecf6270a
[scilab.git] / scilab / modules / umfpack / src / c / taucs_scilab.c
1 /*
2  *   This big file contains all the taucs routines needed for the
3  *   scilab interface onto the taucs supernodal multi-frontal choleski 
4  *   solver :
5  *        - internal structures used by the solver 
6  *        - genmmd interface
7  *        - create / free taucs matrices routines
8  *        - permutations routines
9  *        - symbolic analysis routines
10  *        - super nodal multi frontal factor routines
11  *        - super nodal multi frontal solve routines
12  *        - various utilities
13  *
14  *    All this code is from different files of the TAUCS solver
15  *    lib from Sivan Toledo (see the TAUCS_license.txt file).
16  *
17  *    Very minor modifs have been done (Bruno) :
18  *        - replace the _ for calling fortran routine by the 
19  *          scilab macro F2C  (call to genmmd and calls to blas
20  *          routines)
21  *        - the print messages that I have stayed are redirected
22  *          to the scilab window with sciprint
23  *        - I have added a small routine to get cnz
24  *
25  *    Thanks to Sivan Toledo
26  */
27
28 /*************************************************************/
29 /*                                                           */
30 /*************************************************************/
31
32 #include <stdlib.h>
33 #include <assert.h>
34 #include <math.h>
35 #include "taucs_scilab.h"
36 #include "machine.h"   /* F2C macro */
37 #include "MALLOC.h"
38 #include "BOOL.h"
39 #include "sciprint.h"
40 #include "localization.h"
41
42 #ifdef _MSC_VER
43 #define alloca(x) MALLOC(x)
44 #endif
45
46 #define BLAS_FLOPS_CUTOFF  1000.0
47 #define SOLVE_DENSE_CUTOFF 5
48
49 extern int C2F(genmmd)();
50 extern int C2F(dpotrf)();
51 extern int C2F(dtrsm)();
52 extern int C2F(dsyrk)();
53 extern int C2F(dgemm)();
54 /*************************************************************/
55 /* structures                                                */
56 /*************************************************************/
57
58 typedef struct {
59   int     sn_size;
60   int     n;
61   int*    rowind;
62
63   int     up_size;
64   int*    sn_vertices;
65   int*    up_vertices;
66   double* f1;
67   double* f2;
68   double* u;
69 } supernodal_frontal_matrix;
70
71 typedef struct {
72   char    uplo;     /* 'u' for upper, 'l' for lower, ' ' don't know; prefer lower. */
73   int     n;        /* size of matrix */
74   int     n_sn;     /* number of supernodes */
75
76   int* parent;      /* supernodal elimination tree */
77   int* first_child; 
78   int* next_child;
79
80   int* sn_size;     /* size of supernodes (diagonal block) */
81   int* sn_up_size;  /* size of subdiagonal update blocks   */
82   int** sn_struct;  /* row structure of supernodes         */
83
84   int* sn_blocks_ld;  /* lda of supernode blocks */
85   double** sn_blocks; /* supernode blocks        */
86     
87   int* up_blocks_ld;  /* lda of update blocks    */
88   double** up_blocks; /* update blocks           */
89 } supernodal_factor_matrix;
90
91 /*************************************************************/
92 /* for qsort                                                 */
93 /*************************************************************/
94
95 static int* compare_indirect_map;
96 static int compare_indirect_ints(const void* vx, const void* vy)
97 {
98   int* ix = (int*)vx;
99   int* iy = (int*)vy;
100   if (compare_indirect_map[*ix] < compare_indirect_map[*iy]) return -1;
101   if (compare_indirect_map[*ix] > compare_indirect_map[*iy]) return  1;
102   return 0;
103 }
104
105 /*********************************************************/
106 /* Interface to MMD                                      */
107 /*********************************************************/
108
109 void taucs_ccs_genmmd(taucs_ccs_matrix* m, int** perm, int** invperm)
110 {
111   int  n, maxint, delta, nofsub;
112   int* xadj;
113   int* adjncy;
114   int* invp;
115   int* prm;
116   int* dhead;
117   int* qsize;
118   int* llist;
119   int* marker;
120
121   int* len;
122   int* next;
123
124   int  nnz,i,j,ip;
125   
126   if (!(m->flags & TAUCS_SYMMETRIC)) 
127   {
128         sciprint("%s: %s","taucs_ccs_genmmd",_("GENMMD ordering only works on symmetric matrices.\n"));
129     *perm    = NULL;
130     *invperm = NULL;
131     return;
132   }
133   /* this routine may actually work on UPPER as well */
134   if (!(m->flags & TAUCS_LOWER)) 
135   {
136         sciprint("%s: %s","taucs_ccs_genmmd",_("The lower part of the matrix must be represented.\n"));
137     *perm    = NULL;
138     *invperm = NULL;
139     return;
140   }
141
142   *perm    = NULL;
143   *invperm = NULL;
144
145   n   = m->n;
146   nnz = (m->colptr)[n];
147   
148   /* I copied the value of delta and the size of */
149   /* from SuperLU. Sivan                         */
150
151   delta = 1; /* DELTA is a parameter to allow the choice of nodes
152                 whose degree <= min-degree + DELTA. */
153
154   maxint = 32000;
155
156   assert(sizeof(int) == 4);
157   maxint = 2147483647; /* 2**31-1, for 32-bit only! */
158
159   xadj   = (int*) MALLOC((n+1)     * sizeof(int));
160   adjncy = (int*) MALLOC((2*nnz-n) * sizeof(int));
161   invp   = (int*) MALLOC((n+1)     * sizeof(int));
162   prm    = (int*) MALLOC(n         * sizeof(int));
163   dhead  = (int*) MALLOC((n+1)     * sizeof(int));
164   qsize  = (int*) MALLOC((n+1)     * sizeof(int));
165   llist  = (int*) MALLOC(n         * sizeof(int));
166   marker = (int*) MALLOC(n         * sizeof(int));
167
168   if (!xadj || !adjncy || !invp || !prm 
169       || !dhead || !qsize || !llist || !marker) 
170   {
171     FREE(xadj  );
172     FREE(adjncy);
173     FREE(invp  );
174     FREE(prm   );
175     FREE(dhead );
176     FREE(qsize );
177     FREE(llist );
178     FREE(marker);
179     return;
180   }
181
182   len  = dhead; /* we reuse space */
183   next = qsize; /* we reuse space */
184
185   for (i=0; i<n; i++) len[i] = 0;
186
187   for (j=0; j<n; j++)
188   {
189         for (ip = (m->colptr)[j]; ip < (m->colptr)[j+1]; ip++) 
190     {
191                 i = (m->rowind)[ip];
192                 if (i != j) 
193                 {
194                         len[i] ++;
195                         len[j] ++;
196                 }
197     }
198   }
199   xadj[0] = 1;
200   for (i=1; i<=n; i++) xadj[i] = xadj[i-1] + len[i-1];
201
202   
203   /* use degree as a temporary */
204   for (i=0; i<n; i++) next[i] = xadj[i] - 1;
205
206   for (j=0; j<n; j++)
207   {
208     for (ip = (m->colptr)[j]; ip < (m->colptr)[j+1]; ip++) 
209     {
210                 i = (m->rowind)[ip];
211                 assert( next[i] < 2*nnz-n );
212                 assert( next[j] < 2*nnz-n );
213                 if (i != j) 
214                 {
215                         adjncy[ next[i] ] = j+1;
216                         adjncy[ next[j] ] = i+1;
217                         next[i] ++;
218                         next[j] ++;
219                 }
220     }
221   }
222
223   F2C(genmmd)(&n,
224               xadj, adjncy,
225               invp,prm,
226               &delta,
227               dhead,qsize,llist,marker,
228               &maxint,&nofsub);
229
230   FREE(marker);
231   FREE(llist );
232   FREE(qsize );
233   FREE(dhead );
234   FREE(xadj  );
235   FREE(adjncy);
236
237   for (i=0; i<n; i++) prm[i] --;
238   for (i=0; i<n; i++) invp[ prm[i] ] = i;
239
240   *perm    = prm;
241   *invperm = invp;
242 }
243
244 /*********************************************************/
245 /* CCS : create free taucs matrix routines               */
246 /*********************************************************/
247
248 taucs_ccs_matrix* 
249 taucs_ccs_create(int m, int n, int nnz)
250 {
251   taucs_ccs_matrix* matrix;
252
253   matrix = (taucs_ccs_matrix*) MALLOC(sizeof(taucs_ccs_matrix));
254   if (!matrix) 
255   { 
256         sciprint(_("%s: No more memory.\n"),"taucs_ccs_create");
257     return NULL; 
258   }
259   matrix->flags = 0;
260   matrix->n = n;
261   matrix->m = m;
262   matrix->colptr = (int*)    MALLOC((n+1) * sizeof(int));
263   matrix->rowind = (int*)    MALLOC(nnz   * sizeof(int));
264   matrix->values = (double*) MALLOC(nnz   * sizeof(double));
265   if (!(matrix->colptr) || !(matrix->rowind)) 
266   {
267     sciprint(_("%s: No more memory (n=%d, nnz=%d).\n"),"taucs_ccs_create",n,nnz);
268     FREE(matrix->colptr); FREE(matrix->rowind); FREE(matrix->values);
269     FREE (matrix);
270     return NULL; 
271   }
272
273   return matrix;
274
275
276 void taucs_ccs_free(taucs_ccs_matrix* matrix)
277 {
278   FREE(matrix->rowind);
279   FREE(matrix->colptr);
280   FREE(matrix->values);
281   FREE(matrix);
282 }
283
284 /*********************************************************/
285 /* CCS : permutation routines                            */
286 /*********************************************************/
287
288 taucs_ccs_matrix* 
289 taucs_ccs_permute_symmetrically(taucs_ccs_matrix* A, int* perm, int* invperm)
290 {
291   taucs_ccs_matrix* PAPT;
292   int n;
293   int nnz;
294   /* int* colptr; */
295   int* len;
296   int i,j,ip,I,J;
297   double AIJ;
298
299   /*assert(A->flags & TAUCS_SYMMETRIC);*/
300   assert(A->flags & TAUCS_LOWER);
301
302   n   = A->n;
303   nnz = (A->colptr)[n];
304
305   PAPT = taucs_ccs_create(n,n,nnz);
306   /*PAPT->flags = TAUCS_SYMMETRIC | TAUCS_LOWER;*/
307   PAPT->flags = A->flags;
308
309   len    = (int*) MALLOC(n * sizeof(int));
310   //colptr = (int*) MALLOC(n * sizeof(int));
311
312   for (j=0; j<n; j++) len[j] = 0;
313
314   for (j=0; j<n; j++) 
315   {
316     for (ip = (A->colptr)[j]; ip < (A->colptr)[j+1]; ip++) 
317         {
318       /*i = (A->rowind)[ip] - (A->indshift);*/
319       i = (A->rowind)[ip];
320
321       I = invperm[i];
322       J = invperm[j];
323
324       if (I < J) 
325           {
326                 int T = I; 
327                 J = T;
328       }
329       len[J] ++;
330     }
331   }
332
333   (PAPT->colptr)[0] = 0;
334   for (j=1; j<=n; j++) (PAPT->colptr)[j] = (PAPT->colptr)[j-1] + len[j-1];
335
336   for (j=0; j<n; j++) len[j] = (PAPT->colptr)[j];
337   
338   for (j=0; j<n; j++) 
339   {
340     for (ip = (A->colptr)[j]; ip < (A->colptr)[j+1]; ip++) 
341         {
342       /*i   = (A->rowind)[ip] - (A->indshift);*/
343       i   = (A->rowind)[ip];
344       AIJ = (A->values)[ip];
345
346       I = invperm[i];
347       J = invperm[j];
348
349       if (I < J) 
350           {
351                 int T = I; 
352                 I = J;
353                 J = T;
354       }
355
356       /*(PAPT->rowind)[ len[J] ] = I + (PAPT->indshift);*/
357       (PAPT->rowind)[ len[J] ] = I;
358       (PAPT->values)[ len[J] ] = AIJ;
359
360       len[J] ++;
361     }
362   }
363
364   if (len) {FREE(len); len = NULL;}
365
366   return PAPT;
367 }
368
369 void taucs_vec_permute(int n, double v[], double pv[], int p[])
370 {
371   int i;
372   for (i=0; i<n; i++) pv[i] = v[p[i]];
373
374
375 void taucs_vec_ipermute(int n, double pv[], double v[], int invp[])
376 {
377   int i;
378   for (i=0; i<n; i++) v[invp[i]] = pv[i];
379
380
381
382 /*************************************************************/
383 /* create and free the factor object                         */
384 /*************************************************************/
385
386 static supernodal_factor_matrix*
387 multifrontal_supernodal_create(void)
388 {
389   supernodal_factor_matrix* L;
390   
391   L = (supernodal_factor_matrix*) MALLOC(sizeof(supernodal_factor_matrix));
392   if (!L) return NULL;
393   L->uplo      = 'l';
394   L->n         = -1; /* unused */
395
396   L->sn_struct   = NULL;
397   L->sn_size     = NULL;
398   L->sn_up_size  = NULL;
399   L->parent      = NULL;
400   L->first_child = NULL;
401   L->next_child  = NULL;
402   L->sn_blocks_ld  = NULL;
403   L->sn_blocks     = NULL;
404   L->up_blocks_ld  = NULL;
405   L->up_blocks     = NULL;
406
407   return L;
408 }
409
410 void taucs_supernodal_factor_free(void* vL)
411 {
412   supernodal_factor_matrix* L = (supernodal_factor_matrix*) vL;
413   int sn;
414   
415   FREE(L->parent);
416   FREE(L->first_child);
417   FREE(L->next_child);
418
419   FREE(L->sn_size);
420   FREE(L->sn_up_size);
421   FREE(L->sn_blocks_ld);
422   FREE(L->up_blocks_ld);
423   if (L->sn_struct)   
424     for (sn=0; sn<L->n_sn; sn++)
425       FREE(L->sn_struct[sn]);
426
427   if (L->sn_blocks)   
428     for (sn=0; sn<L->n_sn; sn++)
429       FREE(L->sn_blocks[sn]);
430
431   if (L->up_blocks)   
432     for (sn=0; sn<L->n_sn; sn++)
433       FREE(L->up_blocks[sn]);
434   /*
435   for (sn=0; sn<L->n_sn; sn++) 
436   {
437     FREE(L->sn_struct[sn]);
438     FREE(L->sn_blocks[sn]);
439     FREE(L->up_blocks[sn]);
440         }*/
441
442   FREE(L->sn_struct);
443   FREE(L->sn_blocks);
444   FREE(L->up_blocks);
445
446   FREE(L);
447 }
448
449 /*************************************************************/
450 /* create and free frontal matrices                          */
451 /*************************************************************/
452
453 static supernodal_frontal_matrix* 
454 supernodal_frontal_create(int* firstcol_in_supernode,
455                           int sn_size,
456                           int n, 
457                           int* rowind)
458 {
459   supernodal_frontal_matrix* tmp;
460   int iDontCheckf1 = 0;
461   int iDontCheckf2 = 0;
462   int iDontChecku = 0;
463
464   tmp = (supernodal_frontal_matrix*)MALLOC(sizeof(supernodal_frontal_matrix));
465   if(tmp==NULL) return NULL;
466
467   tmp->sn_size = sn_size;
468   tmp->n = n;
469
470   tmp->rowind = rowind;
471
472   tmp->n = n;
473   tmp->sn_size = sn_size;
474   tmp->up_size = n-sn_size;
475
476   tmp->sn_vertices = rowind;
477   tmp->up_vertices = rowind + sn_size;
478
479   if (tmp->sn_size > 0)
480   {
481           tmp->f1 = (double*)CALLOC((tmp->sn_size)*(tmp->sn_size),sizeof(double));
482   }
483   else
484   {
485         iDontCheckf1 = 1;
486         tmp->f1 = NULL;
487   }
488
489   if (tmp->sn_size  > 0 && tmp->up_size > 0)
490   {
491           tmp->f2 = (double*)CALLOC((tmp->up_size)*(tmp->sn_size),sizeof(double));
492   }
493   else
494   {
495         iDontCheckf2 = 1;
496         tmp->f2 = NULL;
497   }
498   
499   if (tmp->up_size > 0)
500   {
501           tmp->u  = (double*)CALLOC((tmp->up_size)*(tmp->up_size),sizeof(double));
502   }
503   else
504   {
505         iDontChecku = 1;
506         tmp->u = NULL;
507   }
508
509   /*check allocation only if size is > 0 */
510   if( (tmp->f1 == NULL && iDontCheckf1 == 0) || 
511           (tmp->f2 == NULL && iDontCheckf2 == 0) || 
512           (tmp->u  == NULL && iDontChecku  == 0)) 
513   {
514           FREE(tmp->u);
515           FREE(tmp->f1);
516           FREE(tmp->f2);
517           FREE(tmp);
518           return NULL;
519   }
520
521   assert(tmp);
522   return tmp;
523 }
524
525 static void supernodal_frontal_free(supernodal_frontal_matrix* to_del)
526 {
527   /* f1 and f2 are moved to the factor */
528   FREE(to_del->u);
529   FREE(to_del);
530 }
531
532
533 /*************************************************************/
534 /* factor a frontal matrix                                   */
535 /*************************************************************/
536
537 static int
538 multifrontal_supernodal_front_factor(int sn,
539                                      int* firstcol_in_supernode,
540                                      int sn_size,
541                                      taucs_ccs_matrix* A,
542                                      supernodal_frontal_matrix* mtr,
543                                      int* bitmap,
544                                      supernodal_factor_matrix* snL)
545 {
546   int i,j;
547   int* ind;
548   double* re;
549   int INFO = -1;
550   double done      =  1.0;
551   double dminusone = -1.0;
552
553   /* creating transform for real indices */
554   for(i=0;i<mtr->sn_size;i++) bitmap[mtr->sn_vertices[i]] = i;
555   for(i=0;i<mtr->up_size;i++) bitmap[mtr->up_vertices[i]] = mtr->sn_size + i;
556
557   /* adding sn_size column of A to first sn_size column of frontal matrix */
558   for(j=0;j<(mtr->sn_size);j++) 
559   {
560     ind = &(A->rowind[A->colptr[*(firstcol_in_supernode+j)]]);
561     re  = &(A->values[A->colptr[*(firstcol_in_supernode+j)]]); 
562     for(i=0;i < A->colptr[*(firstcol_in_supernode+j)+1]- A->colptr[*(firstcol_in_supernode+j)]; i++) 
563         {
564       if (bitmap[ind[i]] < mtr->sn_size)
565           {
566                 mtr->f1[ (mtr->sn_size)*j + bitmap[ind[i]]] += re[i];
567           }
568       else
569           {
570                 mtr->f2[ (mtr->up_size)*j + bitmap[ind[i]] - mtr->sn_size] += re[i];
571           }
572     }
573   }
574
575   /* we use the BLAS through the Fortran interface */
576
577   /* solving of lower triangular system for L */
578   if (mtr->sn_size)
579   {
580     F2C(dpotrf)("LOWER",
581                 &(mtr->sn_size),
582                 mtr->f1,&(mtr->sn_size),
583                 &INFO);
584   }
585
586   if (INFO && INFO != -1)  /* INFO = -1 if not initialized */ 
587   {
588         sciprint(_("    CC^T Factorization: Matrix is not positive definite.\n"));
589     sciprint(_("                        nonpositive pivot in column %d\n"),mtr->sn_vertices[INFO-1]);
590     return -1;
591   }
592
593   /* getting completion for found columns of L */
594   if (mtr->up_size && mtr->sn_size)
595   {
596     F2C(dtrsm)("Right",
597                "Lower",
598                "Transpose",
599                "No unit diagonal",
600                &(mtr->up_size),&(mtr->sn_size),
601                &done,
602                mtr->f1,&(mtr->sn_size),
603                mtr->f2,&(mtr->up_size));
604   }
605
606   (snL->sn_blocks   )[sn] = mtr->f1;
607   (snL->sn_blocks_ld)[sn] = mtr->sn_size;
608
609   (snL->up_blocks   )[sn] = mtr->f2;
610   (snL->up_blocks_ld)[sn] = mtr->up_size;
611
612   /* computation of updated part of frontal matrix */
613   if (mtr->up_size && mtr->sn_size)
614   {
615     F2C(dsyrk)("Lower",
616                "No Transpose",
617                &(mtr->up_size),&(mtr->sn_size),
618                &dminusone,
619                mtr->f2,&(mtr->up_size),
620                &done,
621                mtr->u, &(mtr->up_size));
622   }
623
624   return 0;
625  }
626
627 /*************************************************************/
628 /* extend-add                                                */
629 /*************************************************************/
630
631 static void 
632 multifrontal_supernodal_front_extend_add(
633                                          supernodal_frontal_matrix* parent_mtr,
634                                          supernodal_frontal_matrix* my_mtr,
635                                          int* bitmap)
636 {
637   int j,i,parent_i,parent_j;
638   double v;
639
640   for(i=0;i<parent_mtr->sn_size;i++) bitmap[parent_mtr->sn_vertices[i]] = i;
641   for(i=0;i<parent_mtr->up_size;i++) bitmap[parent_mtr->up_vertices[i]] = (parent_mtr->sn_size)+i;
642
643   /* extend add operation for update matrix */
644   for(j=0;j<my_mtr->up_size;j++) 
645   {
646     for(i=j;i<my_mtr->up_size;i++) 
647         {
648       parent_j = bitmap[ my_mtr->up_vertices[j] ];
649       parent_i = bitmap[ my_mtr->up_vertices[i] ];
650       /* we could skip this if indices were sorted */
651       if (parent_j>parent_i) 
652           {
653                 int tmp = parent_j;
654                 parent_j = parent_i;
655                 parent_i = tmp;
656       }
657
658       v = (my_mtr->u)[(my_mtr->up_size)*j+i];
659
660       if (parent_j < parent_mtr->sn_size) 
661           {
662                 if (parent_i < parent_mtr->sn_size) 
663                 {
664                         (parent_mtr->f1)[ (parent_mtr->sn_size)*parent_j + parent_i] += v;
665                 } 
666                 else 
667                 {
668                         (parent_mtr->f2)[ (parent_mtr->up_size)*parent_j + (parent_i-parent_mtr->sn_size)] += v;
669                 }
670       } 
671           else 
672           {
673                 (parent_mtr->u)[ (parent_mtr->up_size)*(parent_j-parent_mtr->sn_size) + (parent_i-parent_mtr->sn_size)] += v;
674       }
675     }
676   }
677 }
678
679 /*************************************************************/
680 /* symbolic elimination                                      */
681 /*************************************************************/
682
683 /* UNION FIND ROUTINES */
684
685 static int uf_makeset(int* uf, int i)        { uf[i] = i; return i; }
686 static int uf_find   (int* uf, int i)        
687
688         if (uf[i] != i) uf[i] = uf_find(uf,uf[i]); 
689         return uf[i]; 
690 }
691
692 static int uf_union  (int* uf, int s, int t) {
693   
694   if (uf_find(uf,s) < uf_find(uf,t)) 
695   {
696     uf[uf_find(uf,s)] = uf_find(uf,t); 
697     return (uf_find(uf,t)); 
698   }
699   else
700   {
701     uf[uf_find(uf,s)] = uf_find(uf,t); 
702     return (uf_find(uf,t)); 
703   }
704 }
705
706 static
707 void recursive_postorder(int  j,
708                          int  first_child[],
709                          int  next_child[],
710                          int  postorder[],
711                          int  ipostorder[],
712                          int* next)
713 {
714   int c;
715   for (c=first_child[j]; c != -1; c = next_child[c]) 
716   {
717     recursive_postorder(c,first_child,next_child,
718                         postorder,ipostorder,next);
719   }
720
721   if (postorder)  postorder [*next] = j;
722   if (ipostorder) ipostorder[j] = *next;
723   (*next)++;
724 }
725
726 #define GILBERT_NG_PEYTON_ANALYSIS_SUP
727
728 /* in a few tests the supernodal version seemed slower */
729 #undef GILBERT_NG_PEYTON_ANALYSIS_SUP
730
731 static int ordered_uf_makeset(int* uf, int i)
732
733   uf[i] = i; 
734   return i; 
735 }
736 static int ordered_uf_find   (int* uf, int i) 
737
738   if (uf[i] != i) uf[i] = uf_find(uf,uf[i]); 
739   return uf[i]; 
740 }
741 static int ordered_uf_union  (int* uf, int s, int t) 
742 {
743   assert(uf[t] == t);
744   assert(uf[s] == s);
745   assert(t > s);
746   if (t > s) 
747   {
748     uf[s] = t; 
749     return t; 
750   }
751   else
752   {
753     uf[t] = s;
754   }
755         return s; 
756 }
757
758 static void 
759 tree_level(int j,
760            int isroot, 
761            int first_child[],
762            int next_child[],
763            int level[],
764            int level_j)
765 {
766   int c;
767   if (!isroot) level[j] = level_j;
768   for (c=first_child[j]; c != -1; c = next_child[c]) 
769   {
770     tree_level(c,
771                FALSE,
772                first_child,
773                next_child,
774                level,
775                level_j+1);
776   }
777 }
778
779 static void
780 tree_first_descendant(int j,
781                       int isroot, 
782                       int first_child[],
783                       int next_child[],
784                       int ipostorder[],
785                       int first_descendant[])
786 {
787   int c;
788   int fd = ipostorder[j];
789   for (c=first_child[j]; c != -1; c = next_child[c]) 
790   {
791     tree_first_descendant(c,
792                           FALSE,
793                           first_child,
794                           next_child,
795                           ipostorder,
796                           first_descendant);
797     if (first_descendant[c] < fd) fd = first_descendant[c]; 
798   }
799   if (!isroot) first_descendant[j] = fd;
800 }
801
802 int taucs_ccs_etree(taucs_ccs_matrix* A,
803                 int* parent,
804                 int* l_colcount,
805                 int* l_rowcount,
806                 int* l_nnz)
807 {
808   int* prev_p;
809   int* level;
810   int* l_cc;
811   int* l_rc;
812   int* wt;
813
814   int  n = A->n;
815   int  u,pprime;
816   int  ju;
817   int* postorder;
818   int* ipostorder;
819   int  *first_child,*next_child;
820
821   int i,j,k,ip,jp,kp;
822   int nnz,jnnz;
823   int* uf;
824   int* rowptr;
825   int* colind;
826   int* rowcount;
827   int* realroot;
828
829   /* we need the row structures for the lower triangle */
830
831   nnz = (A->colptr)[n];
832   
833   uf       = (int*)MALLOC(n     * sizeof(int));
834   rowcount = (int*)MALLOC((n+1) * sizeof(int));
835   rowptr   = (int*)MALLOC((n+1) * sizeof(int));
836   colind   = (int*)MALLOC(nnz   * sizeof(int));
837
838   for (i=0; i <=n; i++) rowcount[i] = 0;
839   for (j=0; j < n; j++) {
840     jnnz = (A->colptr)[j+1] - (A->colptr)[j];
841     for (ip=0; ip<jnnz; ip++) {
842       i = (A->rowind)[(A->colptr)[j] + ip];
843       if (j < i) rowcount[i]++;
844     }
845   }
846
847   ip = 0;
848   for (i=0; i <= n; i++) {
849     int next_ip = ip + rowcount[i];
850     rowcount[i] = ip;
851     rowptr  [i] = ip;
852     ip = next_ip;
853   }
854
855   for (j=0; j < n; j++) {
856     jnnz = (A->colptr)[j+1] - (A->colptr)[j];
857     for (ip=0; ip<jnnz; ip++) {
858       i = (A->rowind)[(A->colptr)[j] + ip];
859       if (i==j) continue;
860       assert( rowcount[i] < rowptr[i+1] );
861       colind[ rowcount[i] ] = j;
862       rowcount[i]++;
863     }
864   }
865
866   /* now compute the etree */
867
868   {
869     int u2,t,vroot;
870     realroot = rowcount; /* reuse space */
871
872     for (i=0; i<n; i++) {
873       uf_makeset(uf,i);
874       realroot[i] = i;
875       parent[i] = n;
876       vroot = i;
877       for (kp=rowptr[i]; kp<rowptr[i+1]; kp++) {
878         k=colind[kp];
879         u2 = uf_find(uf,k);
880         t = realroot[u2];
881         if (parent[t] == n && t != i) {
882           parent[t] = i;
883           vroot = uf_union(uf,vroot,u2);
884           realroot[vroot] = i;
885         }
886       }
887     }
888   }
889
890   FREE(colind);
891   FREE(rowptr);
892   FREE(rowcount);
893   
894   /* compute column counts */
895
896   if (l_colcount || l_rowcount || l_nnz) {
897     int* l_nz;
898     int  tmp;
899     int  u2,p,q;
900
901     first_child = (int*)MALLOC((n+1)     * sizeof(int));
902     next_child  = (int*)MALLOC((n+1)     * sizeof(int));
903     postorder   = (int*)MALLOC(n     * sizeof(int));
904     ipostorder  = (int*)MALLOC(n     * sizeof(int));
905     wt  = (int*)MALLOC(n     * sizeof(int));
906     level = (int*)MALLOC(n     * sizeof(int));
907     prev_p  = (int*)MALLOC(n     * sizeof(int));
908
909 #ifdef GILBERT_NG_PEYTON_ANALYSIS_SUP
910     prev_nbr  = MALLOC(n     * sizeof(int));
911     first_descendant  = MALLOC(n     * sizeof(int));
912 #endif
913
914
915     /* compute the postorder */
916     
917     for (j=0; j<=n; j++) first_child[j] = -1;
918     for (j=n-1; j>=0; j--) {
919       next_child[j] = first_child[parent[j]];
920       first_child[parent[j]] = j;
921     }
922     
923     {
924       int next = 0;
925       recursive_postorder(n,first_child,next_child,
926                           postorder,
927                           ipostorder,&next);
928     }
929     
930     /* we allocate scratch vectors to avoid conditionals */
931     /* in the inner loop.                                */
932
933     if (l_colcount) l_cc = l_colcount;
934     else            l_cc = (int*) MALLOC(n*sizeof(int));
935     if (l_rowcount) l_rc = l_rowcount;
936     else            l_rc = (int*) MALLOC(n*sizeof(int));
937     if (l_nnz)      l_nz = l_nnz;
938     else            l_nz = &tmp;
939
940     /* sort by postorder of etree */
941     /* compute level, fst_desc */
942
943     tree_level(n,TRUE,first_child,next_child,
944                level,-1);
945     
946     for (u=0; u < n; u++) prev_p  [u] = -1;
947     for (u=0; u < n; u++) l_rc    [u] =  1;
948     for (u=0; u < n; u++) ordered_uf_makeset(uf,u);
949     for (u=0; u < n; u++) {
950       if (first_child[u] == -1)
951         wt[u] = 1; /* leaves     */
952       else
953         wt[u] =  0; /* nonleaves */
954     }
955
956 #ifdef GILBERT_NG_PEYTON_ANALYSIS_SUP
957     for (u=0; u < n; u++) prev_nbr[u] = -1;
958
959     tree_first_descendant(n,TRUE,
960                           first_child,next_child,ipostorder,
961                           first_descendant);
962 #endif
963
964     FREE(first_child);
965     FREE(next_child);
966
967     for (p=0; p<n; p++) {
968       jp = postorder[p];
969       if (parent[jp] != n) wt[parent[jp]] --;
970       for (ip = (A->colptr)[jp]; ip < (A->colptr)[jp+1]; ip++) {
971         ju = (A->rowind)[ip];
972         if (ju==jp) continue; /* we only want proper neighbors */
973 #ifdef GILBERT_NG_PEYTON_ANALYSIS_SUP
974         if (first_descendant[jp] > prev_nbr[u]) {
975 #else
976         if (1) {
977 #endif
978           wt[jp] ++;
979           pprime = prev_p[ju];
980           if (pprime == -1) 
981             l_rc[ju] += level[jp] - level[ju];
982           else {
983             q = ordered_uf_find(uf,pprime);
984             l_rc[ju] += level[jp] - level[q];
985             wt[q] --;
986           }
987           prev_p[ju] = jp;
988         }
989 #ifdef GILBERT_NG_PEYTON_ANALYSIS_SUP
990         prev_nbr[u] = p;
991 #endif
992       }
993       if (parent[jp] != n) {
994         if (!(ipostorder[parent[jp]] > ipostorder[jp])) {
995           sciprint("jp %d parent %d (ipo_j %d ipo_parent %d\n\r",
996                    jp,parent[jp],ipostorder[jp],ipostorder[parent[jp]]);
997         }
998         assert(ipostorder[parent[jp]] > ipostorder[jp]);
999         ordered_uf_union(uf,jp,parent[jp]);
1000       }
1001     }
1002
1003     *l_nz = 0;
1004     for (u2=0; u2<n; u2++) {
1005       l_cc[u2] = wt[u2];
1006       *l_nz += wt[u2];
1007     }
1008     for (u2=0; u2<n; u2++) {
1009       if (parent[u2] != n) {
1010         l_cc[parent[u2]] += l_cc[u2];
1011         *l_nz += l_cc[u2];
1012       }
1013     }
1014
1015     /* free scrtach vectors                              */
1016
1017     if (!l_colcount) FREE(l_cc);
1018     if (!l_rowcount) FREE(l_rc);
1019
1020     /* free other data structures */
1021
1022     FREE(postorder);
1023     FREE(ipostorder);
1024     FREE(wt);
1025     FREE(level);
1026     FREE(prev_p);
1027     
1028 #ifdef GILBERT_NG_PEYTON_ANALYSIS_SUP
1029     FREE(prev_nbr);
1030     FREE(first_descendant);
1031 #endif
1032   }
1033
1034   FREE(uf);
1035   return 0;
1036 }
1037
1038
1039 int 
1040 taucs_ccs_etree_liu(taucs_ccs_matrix* A,
1041                     int* parent,
1042                     int* l_colcount,
1043                     int* l_rowcount,
1044                     int* l_nnz)
1045 {
1046   int n = A->n;
1047   int i,j,k,ip,kp;
1048   int nnz,jnnz;
1049
1050   int* uf;
1051   int* rowptr;
1052   int* colind;
1053
1054   int* rowcount;
1055   int* marker;
1056   int* realroot;
1057
1058   int* l_cc;
1059   int* l_rc;
1060
1061   /* we need the row structures for the lower triangle */
1062
1063   nnz = (A->colptr)[n];
1064   
1065   uf       = (int*)MALLOC(n     * sizeof(int));
1066   rowcount = (int*)MALLOC((n+1) * sizeof(int));
1067   rowptr   = (int*)MALLOC((n+1) * sizeof(int));
1068   colind   = (int*)MALLOC(nnz   * sizeof(int));
1069
1070   for (i=0; i <=n; i++) rowcount[i] = 0;
1071
1072   for (j=0; j < n; j++) {
1073     
1074     jnnz = (A->colptr)[j+1] - (A->colptr)[j];
1075
1076     for (ip=0; ip<jnnz; ip++) {
1077       i = (A->rowind)[(A->colptr)[j] + ip];
1078       if (j < i) rowcount[i]++;
1079     }
1080   }
1081
1082   ip = 0;
1083   for (i=0; i <= n; i++) {
1084     int next_ip = ip + rowcount[i];
1085     rowcount[i] = ip;
1086     rowptr  [i] = ip;
1087     ip = next_ip;
1088   }
1089
1090   for (j=0; j < n; j++) {
1091     jnnz = (A->colptr)[j+1] - (A->colptr)[j];
1092
1093     for (ip=0; ip<jnnz; ip++) {
1094       i = (A->rowind)[(A->colptr)[j] + ip];
1095       if (i==j) continue;
1096       assert( rowcount[i] < rowptr[i+1] );
1097       colind[ rowcount[i] ] = j;
1098       rowcount[i]++;
1099     }
1100   }
1101
1102   /* now compute the etree */
1103
1104   {
1105     int u,t,vroot;
1106     realroot = rowcount; /* reuse space */
1107
1108     for (i=0; i<n; i++) {
1109       uf_makeset(uf,i);
1110       realroot[i] = i;
1111       parent[i] = n;
1112       vroot = i;
1113       for (kp=rowptr[i]; kp<rowptr[i+1]; kp++) {
1114         k=colind[kp];
1115         u = uf_find(uf,k);
1116         t = realroot[u];
1117         if (parent[t] == n && t != i) {
1118           parent[t] = i;
1119           vroot = uf_union(uf,vroot,u);
1120           realroot[vroot] = i;
1121         }
1122       }
1123     }
1124   }
1125
1126   /* compute column counts */
1127
1128   if (l_colcount || l_rowcount || l_nnz) {
1129     int* l_nz;
1130     int  tmp;
1131
1132     /* we allocate scratch vectors to avoid conditionals */
1133     /* in the inner loop.                                */
1134
1135     if (l_colcount) l_cc = l_colcount;
1136     else            l_cc = (int*) MALLOC(n*sizeof(int));
1137     if (l_rowcount) l_rc = l_rowcount;
1138     else            l_rc = (int*) MALLOC(n*sizeof(int));
1139     if (l_nnz)      l_nz = l_nnz;
1140     else            l_nz = &tmp;
1141
1142     marker = rowcount; /* we reuse the space */
1143     
1144     for (j=0; j < n; j++) l_cc[j] = 1;
1145     *l_nz = n;
1146     
1147     for (i=0; i<n; i++) marker[i] = n; /* clear the array */
1148     
1149     for (i=0; i<n; i++) {
1150       l_rc[i] = 1;
1151       marker[ i ] = i;
1152       for (kp=rowptr[i]; kp<rowptr[i+1]; kp++) {
1153         k=colind[kp];
1154         j=k;
1155         while (marker[j] != i) {
1156           l_cc[j]++;
1157           l_rc[i]++;
1158           (*l_nz)++;
1159           marker[j] = i;
1160           j = parent[j];
1161         }
1162       }
1163     }
1164
1165     /* free scrtach vectors                              */
1166
1167     if (!l_colcount) FREE(l_cc);
1168     if (!l_rowcount) FREE(l_rc);
1169   }
1170
1171   FREE(colind);
1172   FREE(rowptr);
1173   FREE(rowcount);
1174   FREE(uf);
1175   return 0;
1176 }
1177
1178 static void
1179 recursive_symbolic_elimination(int            j,
1180                                taucs_ccs_matrix* A,
1181                                int            first_child[],
1182                                int            next_child[],
1183                                int*           n_sn,
1184                                int            sn_size[],
1185                                int            sn_up_size[],
1186                                int*           sn_rowind[],
1187                                int            sn_first_child[], 
1188                                int            sn_next_child[], 
1189                                int            rowind[],
1190                                int            column_to_sn_map[],
1191                                int            map[],
1192                                int            do_order,
1193                                int            ipostorder[]
1194                                )
1195 {
1196   int  i,ip,c,c_sn;
1197   int  in_previous_sn;
1198   int  nnz;
1199   
1200   for (c=first_child[j]; c != -1; c = next_child[c]) {
1201     recursive_symbolic_elimination(c,A,
1202                                    first_child,next_child,
1203                                    n_sn,
1204                                    sn_size,sn_up_size,sn_rowind,
1205                                    sn_first_child,sn_next_child,
1206                                    rowind, /* temporary */
1207                                    column_to_sn_map,
1208                                    map,
1209                                    do_order,ipostorder
1210                                    );
1211   }
1212
1213   in_previous_sn = 1;
1214   if (j == A->n) 
1215     in_previous_sn = 0; /* this is not a real column */
1216   else if (first_child[j] == -1) 
1217     in_previous_sn = 0; /* this is a leaf */
1218   else if (next_child[first_child[j]] != -1) 
1219     in_previous_sn = 0; /* more than 1 child */
1220   else { 
1221     /* check that the structure is nested */
1222     /* map contains child markers         */
1223
1224     c=first_child[j];
1225     for (ip=(A->colptr)[j]; ip<(A->colptr)[j+1]; ip++) {
1226       i = (A->rowind)[ip];
1227       in_previous_sn = in_previous_sn && (map[i] == c);
1228     }
1229   }
1230
1231   if (in_previous_sn) {
1232     c = first_child[j];
1233     c_sn = column_to_sn_map[c];
1234     column_to_sn_map[j] = c_sn;
1235
1236     /* swap row indices so j is at the end of the */
1237     /* supernode, not in the update indices       */
1238     for (ip=sn_size[c_sn]; ip<sn_up_size[c_sn]; ip++) 
1239       if (sn_rowind[c_sn][ip] == j) break;
1240     assert(ip<sn_up_size[c_sn]);
1241     sn_rowind[c_sn][ip] = sn_rowind[c_sn][sn_size[c_sn]];
1242     sn_rowind[c_sn][sn_size[c_sn]] = j;
1243
1244     /* mark the nonzeros in the map */
1245     for (ip=sn_size[c_sn]; ip<sn_up_size[c_sn]; ip++) 
1246       map[ sn_rowind[c_sn][ip] ] = j;
1247
1248     sn_size   [c_sn]++;
1249
1250     return;
1251   }
1252
1253   /* we are in a new supernode */
1254
1255   if (j < A->n) {
1256     nnz = 1;
1257     rowind[0] = j;
1258     map[j]    = j;
1259     
1260     for (c=first_child[j]; c != -1; c = next_child[c]) {
1261       c_sn = column_to_sn_map[c];
1262       for (ip=sn_size[c_sn]; ip<sn_up_size[c_sn]; ip++) {
1263         i = sn_rowind[c_sn][ip];
1264         if (i > j && map[i] != j) { /* new row index */
1265           map[i] = j;
1266           rowind[nnz] = i;
1267           nnz++;
1268         }
1269       }
1270     }
1271     
1272     for (ip=(A->colptr)[j]; ip<(A->colptr)[j+1]; ip++) {
1273       i = (A->rowind)[ip];
1274       if (map[i] != j) { /* new row index */
1275         map[i] = j;
1276         rowind[nnz] = i;
1277         nnz++;
1278       }
1279     }
1280   }
1281     
1282   /*printf("children of sn %d: ",*n_sn);*/
1283   for (c=first_child[j]; c != -1; c = next_child[c]) {
1284     c_sn = column_to_sn_map[c];
1285     /*printf("%d ",c_sn);*/
1286     if (c==first_child[j])
1287       sn_first_child[*n_sn] = c_sn;
1288     else {
1289       sn_next_child[ c_sn ] = sn_first_child[*n_sn];
1290       sn_first_child[*n_sn] = c_sn;
1291     }
1292   }
1293   /*printf("\n");*/
1294
1295   if (j < A->n) {
1296     column_to_sn_map[j] = *n_sn;
1297     sn_size   [*n_sn] = 1;
1298     sn_up_size[*n_sn] = nnz;
1299     sn_rowind [*n_sn] = (int*) MALLOC(nnz * sizeof(int));
1300     for (ip=0; ip<nnz; ip++) sn_rowind[*n_sn][ip] = rowind[ip];
1301     if (do_order) {
1302       /* Sivan and Vladimir: we think that we can sort in */
1303       /* column order, not only in etree postorder.       */
1304       /*
1305         radix_sort(sn_rowind [*n_sn],nnz);
1306         qsort(sn_rowind [*n_sn],nnz,sizeof(int),compare_ints);
1307       */
1308       compare_indirect_map = ipostorder;
1309       qsort(sn_rowind [*n_sn],nnz,sizeof(int),compare_indirect_ints);
1310     }
1311     assert(sn_rowind [*n_sn][0] == j);
1312     (*n_sn)++;
1313   }
1314
1315   return;
1316 }
1317
1318 /* count zeros and nonzeros in a supernode to compute the */
1319 /* utility of merging fundamental supernodes.             */
1320
1321 typedef struct {
1322   double zeros;
1323   double nonzeros;
1324 } znz;
1325
1326 static znz
1327 recursive_amalgamate_supernodes(int           sn,
1328                                 int*           n_sn,
1329                                 int            sn_size[],
1330                                 int            sn_up_size[],
1331                                 int*           sn_rowind[],
1332                                 int            sn_first_child[], 
1333                                 int            sn_next_child[], 
1334                                 int            rowind[],
1335                                 int            column_to_sn_map[],
1336                                 int            map[],
1337                                 int            do_order,
1338                                 int            ipostorder[]
1339                                 )
1340 {
1341   int  i,ip,c_sn,gc_sn;
1342   int  nnz;
1343   int  nchildren; /* number of children, child index */
1344   znz* c_znz;
1345   znz  sn_znz, merged_znz;
1346
1347   int new_sn_size, new_sn_up_size;
1348
1349   sn_znz.zeros    = 0.0;
1350   sn_znz.nonzeros = (double) (((sn_up_size[sn] - sn_size[sn]) * sn_size[sn]) 
1351                               + (sn_size[sn] * (sn_size[sn] + 1))/2);
1352
1353   if (sn_first_child[sn] == -1) { /* leaf */
1354     return sn_znz;
1355   }
1356
1357   nchildren = 0;
1358   for (c_sn=sn_first_child[sn]; c_sn != -1; c_sn = sn_next_child[c_sn])
1359     nchildren++;
1360
1361   c_znz = (znz*) alloca(nchildren * sizeof(znz));
1362
1363   //printf("supernode %d out of %d\n",sn,*n_sn);
1364
1365   /* merge the supernode with its children! */
1366
1367   i = 0;
1368   for (c_sn=sn_first_child[sn]; c_sn != -1; c_sn = sn_next_child[c_sn]) {
1369     c_znz[i] = 
1370       recursive_amalgamate_supernodes(c_sn,
1371                                       n_sn,
1372                                       sn_size,sn_up_size,sn_rowind,
1373                                       sn_first_child,sn_next_child,
1374                                       rowind, /* temporary */
1375                                       column_to_sn_map,
1376                                       map,
1377                                       do_order,ipostorder
1378                                       );
1379     assert(c_znz[i].zeros + c_znz[i].nonzeros ==
1380            (double) (((sn_up_size[c_sn] - sn_size[c_sn]) * sn_size[c_sn]) 
1381                      + (sn_size[c_sn] * (sn_size[c_sn] + 1))/2 ));
1382     i++;
1383   }
1384
1385   merged_znz.nonzeros = sn_znz.nonzeros;
1386   merged_znz.zeros    = sn_znz.zeros;
1387                    
1388   for (i=0; i<nchildren; i++) {
1389     merged_znz.nonzeros += (c_znz[i]).nonzeros;
1390     merged_znz.zeros    += (c_znz[i]).zeros;
1391   }
1392
1393   /*  printf("supernode %d out of %d (continuing)\n",sn,*n_sn);*/
1394
1395   /* should we merge the supernode with its children? */
1396
1397   nnz = 0;
1398   for (c_sn=sn_first_child[sn]; c_sn != -1; c_sn = sn_next_child[c_sn]) {
1399     for (ip=0; ip<sn_size[c_sn]; ip++) {
1400       i = sn_rowind[c_sn][ip];
1401       assert( map[i] != sn );
1402       map[i] = sn;
1403       rowind[nnz] = i;
1404       nnz++;
1405     }
1406   }
1407
1408   for (ip=0; ip<sn_size[sn]; ip++) {
1409     i = sn_rowind[sn][ip];
1410     assert( map[i] != sn );
1411     map[i] = sn;
1412     rowind[nnz] = i;
1413     nnz++;
1414   }
1415
1416   new_sn_size = nnz;
1417
1418   for (c_sn=sn_first_child[sn]; c_sn != -1; c_sn = sn_next_child[c_sn]) {
1419     for (ip=sn_size[c_sn]; ip<sn_up_size[c_sn]; ip++) {
1420       i = sn_rowind[c_sn][ip];
1421       if (map[i] != sn) { /* new row index */
1422         map[i] = sn;
1423         rowind[nnz] = i;
1424         nnz++;
1425       }
1426     }
1427   }
1428
1429   for (ip=sn_size[sn]; ip<sn_up_size[sn]; ip++) {
1430     i = sn_rowind[sn][ip];
1431     if (map[i] != sn) { /* new row index */
1432       map[i] = sn;
1433       rowind[nnz] = i;
1434       nnz++;
1435     }
1436   }
1437   
1438   new_sn_up_size = nnz;
1439
1440   if (do_order) {
1441     compare_indirect_map = ipostorder;
1442     qsort(rowind,nnz,sizeof(int),compare_indirect_ints);
1443   }
1444
1445   /* determine whether we should merge the supernode and its children */
1446
1447   {
1448     int n;
1449     double* zcount;
1450
1451     n = 0;
1452     for (ip=0; ip<nnz; ip++) {
1453       i = rowind[ip];
1454       if (i >= n) n = i+1;
1455     }
1456
1457     zcount = (double*) malloc(n * sizeof(double));
1458     assert(zcount);
1459     
1460     for (ip=0; ip<new_sn_size; ip++) {
1461       i = rowind[ip]; assert(i<n);
1462       zcount[i] = (double) (ip+1);
1463     }
1464     for (ip=new_sn_size; ip<new_sn_up_size; ip++) {
1465       i = rowind[ip]; assert(i<n);
1466       zcount[i] = (double) new_sn_size;
1467     }
1468
1469     /*
1470     for (ip=0; ip<new_sn_up_size; ip++) 
1471       printf("row %d zcount = %.0f\n",rowind[ip],zcount[rowind[ip]]);
1472     */
1473     
1474     for (c_sn=sn_first_child[sn]; c_sn != -1; c_sn = sn_next_child[c_sn]) {
1475       for (ip=0; ip<sn_size[c_sn]; ip++) {
1476         i = sn_rowind[c_sn][ip]; assert(i<n);
1477         zcount[i] -= (double) (ip+1);
1478       }
1479       for (ip=sn_size[c_sn]; ip<sn_up_size[c_sn]; ip++) {
1480         i = sn_rowind[c_sn][ip]; assert(i<n);
1481         zcount[i] -= (double) sn_size[c_sn];
1482       }
1483     }
1484
1485     for (ip=0; ip<sn_size[sn]; ip++) {
1486       i = sn_rowind[sn][ip]; assert(i<n);
1487       zcount[i] -= (double) (ip+1);
1488     }
1489     for (ip=sn_size[sn]; ip<sn_up_size[sn]; ip++) {
1490       i = sn_rowind[sn][ip]; assert(i<n);
1491       zcount[i] -= (double) sn_size[sn];
1492     }
1493
1494     /*
1495     for (ip=0; ip<new_sn_up_size; ip++) 
1496       printf("ROW %d zcount = %.0f\n",rowind[ip],zcount[rowind[ip]]);
1497     printf("zeros before merging %.0f\n",merged_znz.zeros);
1498     */
1499     
1500     for (ip=0; ip<new_sn_up_size; ip++) {
1501       i = rowind[ip]; assert(i<n);
1502       assert(zcount[i] >= 0.0);
1503       merged_znz.zeros += zcount[i];
1504     }
1505     FREE( zcount );
1506     /*printf("zeros after merging %.0f\n",merged_znz.zeros);*/
1507
1508     /* voodoo constants (need some kind of a utility function */
1509     if ((new_sn_size < 16)
1510         ||
1511         ((sn_size[sn] < 50) && (merged_znz.zeros < 0.5 * merged_znz.nonzeros))
1512         ||
1513         ((sn_size[sn] < 250) && (merged_znz.zeros < 0.25 * merged_znz.nonzeros))
1514         ||
1515         ((sn_size[sn] < 500) && (merged_znz.zeros < 0.10 * merged_znz.nonzeros))
1516         ||
1517         (merged_znz.zeros < 0.05 * merged_znz.nonzeros)
1518         ) {
1519       /*
1520       taucs_printf("merging sn %d, zeros (%f) vs nonzeros (%f)\n",
1521                    sn,merged_znz.zeros,merged_znz.nonzeros);
1522       */
1523     } else {
1524       /*
1525       taucs_printf("sn %d, too many zeros (%f) vs nonzeros (%f)\n",
1526                    sn,merged_znz.zeros,merged_znz.nonzeros);
1527       printf("returning without merging\n");
1528       */
1529       return sn_znz;
1530     }
1531   }
1532
1533   /* now merge the children lists */
1534
1535   sn_size[sn]    = new_sn_size;
1536   sn_up_size[sn] = new_sn_up_size;
1537   sn_rowind[sn]  = (int*) REALLOC(sn_rowind[sn], 
1538                                   new_sn_up_size * sizeof(int));
1539   for (ip=0; ip<new_sn_up_size; ip++) sn_rowind[sn][ip] = rowind[ip];
1540
1541   //  printf("supernode %d out of %d (merging)\n",sn,*n_sn);
1542
1543   nchildren = 0;
1544   for (c_sn=sn_first_child[sn]; c_sn != -1; c_sn = sn_next_child[c_sn]) {
1545     for (ip=0; ip<sn_size[c_sn]; ip++) {
1546       i = (sn_rowind[c_sn])[ip];
1547       assert(column_to_sn_map[i] == c_sn);
1548       column_to_sn_map[i] = sn;
1549     }
1550
1551     for (gc_sn=sn_first_child[c_sn]; gc_sn != -1; gc_sn = sn_next_child[gc_sn]) {
1552       rowind[nchildren] = gc_sn;
1553       nchildren++;
1554     }
1555   }
1556
1557   /* free the children's rowind vectors */
1558   for (c_sn=sn_first_child[sn]; c_sn != -1; c_sn = sn_next_child[c_sn]) {
1559     FREE( sn_rowind[c_sn] );
1560     sn_rowind[c_sn]  = NULL;
1561     sn_size[c_sn]    = 0;
1562     sn_up_size[c_sn] = 0;
1563   }
1564
1565   sn_first_child[sn] = -1;
1566   for (i=0; i<nchildren; i++) {
1567     sn_next_child[ rowind[i] ] = sn_first_child[sn];
1568     sn_first_child[sn] = rowind[i];
1569   }    
1570
1571   /*
1572   printf("supernode %d out of %d (done)\n",sn,*n_sn);
1573   printf("returning, merging\n");
1574   */
1575   return merged_znz;
1576 }
1577
1578 int      
1579 taucs_ccs_symbolic_elimination(taucs_ccs_matrix* A,
1580                                void* vL,
1581                                int do_order
1582                                )
1583 {
1584   supernodal_factor_matrix* L = (supernodal_factor_matrix*) vL;
1585   int* first_child;
1586   int* next_child;
1587   int j;
1588   int* column_to_sn_map;
1589   int* map;
1590   int* rowind;
1591   int* parent;
1592   int* ipostorder;
1593
1594   L->n         = A->n;
1595   L->sn_struct = (int**)MALLOC((A->n  )*sizeof(int*));
1596   L->sn_size   = (int*) MALLOC((A->n+1)*sizeof(int));
1597   
1598   L->sn_up_size   = (int*) MALLOC((A->n+1)*sizeof(int));
1599   L->first_child = (int*) MALLOC((A->n+1)*sizeof(int));
1600   L->next_child  = (int*) MALLOC((A->n+1)*sizeof(int));
1601   
1602   column_to_sn_map = (int*)MALLOC((A->n+1)*sizeof(int));
1603   map              = (int*) MALLOC((A->n+1)*sizeof(int));
1604
1605   first_child = (int*) MALLOC(((A->n)+1)*sizeof(int));
1606   next_child  = (int*) MALLOC(((A->n)+1)*sizeof(int));
1607     
1608   rowind      = (int*) MALLOC((A->n)*sizeof(int));
1609
1610   /* compute the vertex elimination tree */
1611   parent      = (int*)MALLOC((A->n+1)*sizeof(int));
1612
1613   taucs_ccs_etree(A,parent,NULL,NULL,NULL);
1614
1615   if (0) {
1616     int *cc1,*cc2,*rc1,*rc2;
1617     int *p1;
1618     int nnz1,nnz2;
1619
1620     cc1=(int*)MALLOC((A->n)*sizeof(int));
1621     cc2=(int*)MALLOC((A->n)*sizeof(int));
1622     rc1=(int*)MALLOC((A->n)*sizeof(int));
1623     rc2=(int*)MALLOC((A->n)*sizeof(int));
1624     p1 =(int*)MALLOC((A->n)*sizeof(int));
1625
1626     taucs_ccs_etree_liu(A,parent,cc1,rc1,&nnz1);
1627
1628     taucs_ccs_etree(A,p1,cc2,rc2,&nnz2);
1629
1630     for (j=0; j<(A->n); j++) assert(parent[j]==p1[j]);
1631     for (j=0; j<(A->n); j++) {
1632       if (cc1[j]!=cc2[j]) sciprint("j=%d cc1=%d cc2=%d\n",j,cc1[j],cc2[j]);
1633       assert(cc1[j]==cc2[j]);
1634     }
1635
1636     for (j=0; j<(A->n); j++) {
1637       if (rc1[j]!=rc2[j]) sciprint("j=%d rc1=%d rc2=%d\n",j,rc1[j],rc2[j]);
1638       assert(rc1[j]==rc2[j]);
1639     }
1640
1641     if (nnz1!=nnz2) sciprint("nnz1=%d nnz2=%d\n",nnz1,nnz2);
1642     
1643     FREE(cc1); FREE(cc2); FREE(rc1); FREE(rc2);
1644   }
1645
1646   for (j=0; j <= (A->n); j++) first_child[j] = -1;
1647   for (j = (A->n)-1; j >= 0; j--) {
1648     int p = parent[j];
1649     next_child[j] = first_child[p];
1650     first_child[p] = j;
1651   }
1652   FREE(parent);
1653
1654   ipostorder = (int*)MALLOC((A->n+1)*sizeof(int));
1655   { 
1656     int next = 0;
1657     /*int* postorder = (int*)MALLOC((A->n+1)*sizeof(int));*/
1658     recursive_postorder(A->n,first_child,next_child,
1659                         NULL,
1660                         ipostorder,&next);
1661     /*
1662     printf("ipostorder ");
1663     for (j=0; j <= (A->n); j++) printf("%d ",ipostorder[j]);
1664     printf("\n");
1665     printf(" postorder ");
1666     for (j=0; j <= (A->n); j++) printf("%d ",postorder[j]);
1667     printf("\n");
1668     */
1669   }
1670
1671   L->n_sn = 0;
1672   for (j=0; j < (A->n); j++) map[j] = -1;
1673   for (j=0; j <= (A->n); j++) (L->first_child)[j] = (L->next_child)[j] = -1;
1674   
1675   recursive_symbolic_elimination(A->n,
1676                                  A,
1677                                  first_child,next_child,
1678                                  &(L->n_sn),
1679                                  L->sn_size,L->sn_up_size,L->sn_struct,
1680                                  L->first_child,L->next_child,
1681                                  rowind,
1682                                  column_to_sn_map,
1683                                  map,
1684                                  do_order,ipostorder
1685                                  );
1686
1687   {
1688     double nnz   = 0.0;
1689     double flops = 0.0;
1690     int sn,i,colnnz;
1691     for (sn=0; sn<(L->n_sn); sn++) {
1692       for (i=0, colnnz = (L->sn_up_size)[sn]; 
1693            i<(L->sn_size)[sn]; 
1694            i++, colnnz--) {
1695         flops += ((double)(colnnz) - 1.0) * ((double)(colnnz) + 2.0) / 2.0;
1696         nnz   += (double) (colnnz);
1697       }
1698     }
1699     /*sciprint("\t\tSymbolic Analysis of CC^T: %.2e nonzeros, %.2e flops\n",
1700       nnz, flops); */
1701   }
1702
1703   for (j=0; j < (A->n); j++) map[j] = -1;
1704
1705   (void) recursive_amalgamate_supernodes((L->n_sn) - 1,
1706                                          &(L->n_sn),
1707                                          L->sn_size,L->sn_up_size,L->sn_struct,
1708                                          L->first_child,L->next_child,
1709                                          rowind,
1710                                          column_to_sn_map,
1711                                          map,
1712                                          do_order,ipostorder
1713                                          );
1714
1715
1716   {
1717     double nnz   = 0.0;
1718     double flops = 0.0;
1719     int sn,i,colnnz;
1720     for (sn=0; sn<(L->n_sn); sn++) {
1721       for (i=0, colnnz = (L->sn_up_size)[sn]; 
1722            i<(L->sn_size)[sn]; 
1723            i++, colnnz--) {
1724         flops += ((double)(colnnz) - 1.0) * ((double)(colnnz) + 2.0) / 2.0;
1725         nnz   += (double) (colnnz);
1726       }
1727     }
1728     /*sciprint("\t\tRelaxed  Analysis of LL^T: %.2e nonzeros, %.2e flops\n",
1729       nnz, flops); */
1730   }
1731
1732   /*
1733   {
1734     int i;
1735     printf("c2sn: ");
1736     for (i=0; i<A->n; i++) printf("%d ",column_to_sn_map[i]);
1737     printf("\n");
1738   }
1739   */
1740   
1741
1742
1743   L->sn_blocks_ld  = (int*)MALLOC((L->n_sn) * sizeof(int));
1744   L->sn_blocks     = (double**)CALLOC((L->n_sn), sizeof(double*)); /* so we can free before allocation */
1745   
1746   L->up_blocks_ld  = (int*)MALLOC((L->n_sn) * sizeof(int));
1747   L->up_blocks     = (double**)CALLOC((L->n_sn), sizeof(double*));
1748
1749   FREE(rowind);
1750   FREE(map);
1751   FREE(column_to_sn_map);
1752   FREE(next_child);
1753   FREE(first_child);
1754   FREE(ipostorder);
1755
1756   return 0;
1757 }
1758
1759
1760 /*************************************************************/
1761 /* factor routines                                           */
1762 /*************************************************************/
1763
1764 static supernodal_frontal_matrix*
1765 recursive_multifrontal_supernodal_factor_llt(int sn,       /* this supernode */
1766                                              int is_root,  /* is v the root? */
1767                                              int* bitmap,
1768                                              taucs_ccs_matrix* A,
1769                                              supernodal_factor_matrix* snL,
1770                                              int* fail)
1771 {
1772   supernodal_frontal_matrix* my_matrix=NULL;
1773   supernodal_frontal_matrix* child_matrix=NULL;
1774   int child;
1775   int* v; 
1776   int  sn_size;
1777   int* first_child   = snL->first_child;
1778   int* next_child    = snL->next_child;
1779
1780   v = &( snL->sn_struct[sn][0] );
1781
1782   if (!is_root)
1783     sn_size = snL->sn_size[sn];
1784   else
1785     sn_size = -1;
1786
1787   for (child = first_child[sn]; child != -1; child = next_child[child]) {
1788     child_matrix = 
1789       recursive_multifrontal_supernodal_factor_llt(child,
1790                                                    FALSE,
1791                                                    bitmap,
1792                                                    A,snL,fail);
1793     if (*fail) { 
1794       if (child_matrix) supernodal_frontal_free(child_matrix);
1795       return NULL;
1796     }
1797
1798     if (!is_root) {
1799       if (!my_matrix) {
1800         my_matrix =  supernodal_frontal_create(v,sn_size,
1801                                                snL->sn_up_size[sn],
1802                                                snL->sn_struct[sn]);
1803         if (!my_matrix) {
1804           *fail = TRUE;
1805           supernodal_frontal_free(child_matrix);
1806           return NULL;
1807         }
1808       }
1809
1810       multifrontal_supernodal_front_extend_add(my_matrix,child_matrix,bitmap);
1811       supernodal_frontal_free(child_matrix);
1812     }
1813   }
1814
1815   /* in case we have no children, we allocate now */
1816   if (!is_root && !my_matrix) 
1817     my_matrix =  supernodal_frontal_create(v,sn_size,
1818                                            snL->sn_up_size[sn],
1819                                            snL->sn_struct[sn]);
1820
1821   if(!is_root) {
1822     if (multifrontal_supernodal_front_factor(sn,
1823                                              v,sn_size,
1824                                              A,
1825                                              my_matrix,
1826                                              bitmap,
1827                                              snL)) { 
1828       /* nonpositive pivot */
1829       *fail = TRUE;
1830       supernodal_frontal_free(my_matrix);
1831       return NULL;
1832     }
1833   }
1834   return my_matrix;
1835 }
1836
1837 void* taucs_ccs_factor_llt_mf(taucs_ccs_matrix* A)
1838 {
1839   supernodal_factor_matrix* L;
1840   int* map;
1841   int fail;
1842
1843   L = multifrontal_supernodal_create();
1844
1845   taucs_ccs_symbolic_elimination(A,L,
1846                                  FALSE /* don't sort row indices */ );
1847
1848   map = (int*)MALLOC((A->n+1)*sizeof(int));
1849
1850   fail = FALSE;
1851   recursive_multifrontal_supernodal_factor_llt((L->n_sn),  
1852                                                TRUE, 
1853                                                map,
1854                                                A,L,&fail);
1855
1856   FREE(map);
1857
1858   if (fail) {
1859     taucs_supernodal_factor_free(L);
1860     return NULL;
1861   }
1862
1863   return (void*) L;
1864 }
1865
1866
1867 /*************************************************************/
1868 /* supernodal solve routines                                 */
1869 /*************************************************************/
1870
1871 static void 
1872 recursive_supernodal_solve_l(int sn,       /* this supernode */
1873                              int is_root,  /* is v the root? */
1874                              int* first_child, int* next_child,
1875                              int** sn_struct, int* sn_sizes, int* sn_up_sizes,
1876                              int* sn_blocks_ld,double* sn_blocks[],
1877                              int* up_blocks_ld,double* up_blocks[],
1878                              double x[], double b[],
1879                              double t[])
1880 {
1881   int child;
1882   int  sn_size; /* number of rows/columns in the supernode    */
1883   int  up_size; /* number of rows that this supernode updates */
1884   int    ione = 1;
1885   double done = 1.0;
1886   double dzero = 0.0;
1887   double* xdense;
1888   double* bdense;
1889   double  flops;
1890   int i,j,ip,jp;
1891
1892   for (child = first_child[sn]; child != -1; child = next_child[child]) {
1893     recursive_supernodal_solve_l(child,
1894                                  FALSE,
1895                                  first_child,next_child,
1896                                  sn_struct,sn_sizes,sn_up_sizes,
1897                                  sn_blocks_ld,sn_blocks,
1898                                  up_blocks_ld,up_blocks,
1899                                  x,b,t);
1900   }
1901
1902   if(!is_root) {
1903
1904     sn_size = sn_sizes[sn];
1905     up_size = sn_up_sizes[sn] - sn_sizes[sn];
1906
1907     flops = ((double)sn_size)*((double)sn_size) 
1908       + 2.0*((double)sn_size)*((double)up_size);
1909
1910     if (flops > BLAS_FLOPS_CUTOFF) {
1911       xdense = t;
1912       bdense = t + sn_size;
1913       
1914       for (i=0; i<sn_size; i++)
1915         xdense[i] = b[ sn_struct[ sn ][ i ] ];
1916       for (i=0; i<up_size; i++)
1917         bdense[i] = 0.0;
1918       
1919       F2C(dtrsm)("Left",
1920                  "Lower",
1921                  "No Transpose",
1922                  "No unit diagonal",
1923                  &sn_size,&ione,
1924                  &done,
1925                  sn_blocks[sn],&(sn_blocks_ld[sn]),
1926                  xdense       ,&sn_size);
1927       
1928       if ((up_size > 0) & (sn_size > 0)) {
1929         F2C(dgemm)("No Transpose","No Transpose",
1930                    &up_size, &ione, &sn_size,
1931                    &done,
1932                    up_blocks[sn],&(up_blocks_ld[sn]),
1933                    xdense       ,&sn_size,
1934                    &dzero,
1935                    bdense       ,&up_size);
1936       }
1937       
1938       for (i=0; i<sn_size; i++)
1939         x[ sn_struct[ sn][ i ] ]  = xdense[i];
1940       for (i=0; i<up_size; i++)
1941         b[ sn_struct[ sn ][ sn_size + i ] ] -= bdense[i];
1942
1943     } else if (sn_size > SOLVE_DENSE_CUTOFF) {
1944
1945       xdense = t;
1946       bdense = t + sn_size;
1947       
1948       for (i=0; i<sn_size; i++)
1949         xdense[i] = b[ sn_struct[ sn ][ i ] ];
1950       for (i=0; i<up_size; i++)
1951         bdense[i] = 0.0;
1952       
1953       for (jp=0; jp<sn_size; jp++) {
1954         xdense[jp] = xdense[jp] / sn_blocks[sn][ sn_blocks_ld[sn]*jp + jp];
1955
1956         for (ip=jp+1; ip<sn_size; ip++) {
1957           xdense[ip] -= xdense[jp] * sn_blocks[sn][ sn_blocks_ld[sn]*jp + ip];
1958         }
1959       }
1960
1961       for (jp=0; jp<sn_size; jp++) {
1962         for (ip=0; ip<up_size; ip++) {
1963           bdense[ip] += xdense[jp] * up_blocks[sn][ up_blocks_ld[sn]*jp + ip];
1964         }
1965       }
1966
1967       for (i=0; i<sn_size; i++)
1968         x[ sn_struct[ sn][ i ] ]  = xdense[i];
1969       for (i=0; i<up_size; i++)
1970         b[ sn_struct[ sn ][ sn_size + i ] ] -= bdense[i];
1971       
1972     } else {
1973
1974       for (jp=0; jp<sn_size; jp++) {
1975         j = sn_struct[sn][jp];
1976         x[j] = b[j] / sn_blocks[sn][ sn_blocks_ld[sn]*jp + jp];
1977         for (ip=jp+1; ip<sn_size; ip++) {
1978           i = sn_struct[sn][ip];
1979           b[i] -= x[j] * sn_blocks[sn][ sn_blocks_ld[sn]*jp + ip];
1980         }
1981
1982         for (ip=0; ip<up_size; ip++) {
1983           i = sn_struct[sn][sn_size + ip];
1984           b[i] -= x[j] * up_blocks[sn][ up_blocks_ld[sn]*jp + ip];
1985         }
1986       }
1987
1988     }
1989   }
1990 }
1991
1992 static void 
1993 recursive_supernodal_solve_lt(int sn,       /* this supernode */
1994                               int is_root,  /* is v the root? */
1995                               int* first_child, int* next_child,
1996                               int** sn_struct, int* sn_sizes, int* sn_up_sizes,
1997                               int* sn_blocks_ld,double* sn_blocks[],
1998                               int* up_blocks_ld,double* up_blocks[],
1999                               double x[], double b[],
2000                               double t[])
2001 {
2002   int child;
2003   int  sn_size; /* number of rows/columns in the supernode    */
2004   int  up_size; /* number of rows that this supernode updates */
2005   int    ione = 1;
2006   double done = 1.0;
2007   double dminusone = -1.0;
2008   double* xdense;
2009   double* bdense;
2010   double  flops;
2011   int i,j,ip,jp;
2012
2013   if(!is_root) {
2014
2015     sn_size = sn_sizes[sn];
2016     up_size = sn_up_sizes[sn]-sn_sizes[sn];
2017     
2018     flops = ((double)sn_size)*((double)sn_size) 
2019       + 2.0*((double)sn_size)*((double)up_size);
2020
2021     if (flops > BLAS_FLOPS_CUTOFF) {
2022
2023       bdense = t;
2024       xdense = t + sn_size;
2025       
2026       for (i=0; i<sn_size; i++)
2027         bdense[i] = b[ sn_struct[ sn][ i ] ];
2028       for (i=0; i<up_size; i++)
2029         xdense[i] = x[ sn_struct[sn][sn_size+i] ];
2030       
2031       if ((up_size > 0) & (sn_size > 0))
2032         F2C(dgemm)("Transpose","No Transpose",
2033                    &sn_size, &ione, &up_size,
2034                    &dminusone,
2035                    up_blocks[sn],&(up_blocks_ld[sn]),
2036                    xdense       ,&up_size,
2037                    &done,
2038                    bdense       ,&sn_size);
2039       
2040       F2C(dtrsm)("Left",
2041                  "Lower",
2042                  "Transpose",
2043                  "No unit diagonal",
2044                  &sn_size,&ione,
2045                  &done,
2046                  sn_blocks[sn],&(sn_blocks_ld[sn]),
2047                  bdense       ,&sn_size);
2048       
2049       for (i=0; i<sn_size; i++)
2050         x[ sn_struct[ sn][ i ] ]  = bdense[i];
2051     
2052     } else if (sn_size > SOLVE_DENSE_CUTOFF) {
2053
2054       bdense = t;
2055       xdense = t + sn_size;
2056       
2057       for (i=0; i<sn_size; i++)
2058         bdense[i] = b[ sn_struct[ sn][ i ] ];
2059       for (i=0; i<up_size; i++)
2060         xdense[i] = x[ sn_struct[sn][sn_size+i] ];
2061       
2062       for (ip=sn_size-1; ip>=0; ip--) {
2063         for (jp=0; jp<up_size; jp++) {
2064           bdense[ip] -= xdense[jp] * up_blocks[sn][ up_blocks_ld[sn]*ip + jp];
2065         }
2066       }
2067
2068       for (ip=sn_size-1; ip>=0; ip--) {
2069         for (jp=sn_size-1; jp>ip; jp--) {
2070           bdense[ip] -= bdense[jp] * sn_blocks[sn][ sn_blocks_ld[sn]*ip + jp];
2071         }
2072         bdense[ip] = bdense[ip] / sn_blocks[sn][ sn_blocks_ld[sn]*ip + ip];
2073       }
2074
2075       for (i=0; i<sn_size; i++)
2076         x[ sn_struct[ sn][ i ] ]  = bdense[i];
2077     
2078     } else {
2079
2080       for (ip=sn_size-1; ip>=0; ip--) {
2081         i = sn_struct[sn][ip];
2082
2083         for (jp=0; jp<up_size; jp++) {
2084           j = sn_struct[sn][sn_size + jp];
2085           b[i] -= x[j] * up_blocks[sn][ up_blocks_ld[sn]*ip + jp];
2086         }
2087
2088         for (jp=sn_size-1; jp>ip; jp--) {
2089           j = sn_struct[sn][jp];
2090           b[i] -= x[j] * sn_blocks[sn][ sn_blocks_ld[sn]*ip + jp];
2091         }
2092         x[i] = b[i] / sn_blocks[sn][ sn_blocks_ld[sn]*ip + ip];
2093
2094       }
2095
2096     }
2097   }
2098
2099   for (child = first_child[sn]; child != -1; child = next_child[child]) {
2100     recursive_supernodal_solve_lt(child,
2101                                   FALSE,
2102                                   first_child,next_child,
2103                                   sn_struct,sn_sizes,sn_up_sizes,
2104                                   sn_blocks_ld,sn_blocks,
2105                                   up_blocks_ld,up_blocks,
2106                                   x,b,t);
2107   }
2108 }
2109
2110 int taucs_supernodal_solve_llt(void* vL,
2111                                double* x, double* b)
2112 {
2113   supernodal_factor_matrix* L = (supernodal_factor_matrix*) vL;
2114   double* y;
2115   double* t; /* temporary vector */
2116   int     i;
2117   
2118   y = MALLOC((L->n) * sizeof(double));
2119   t = MALLOC((L->n) * sizeof(double));
2120   if (!y || !t) {
2121     FREE(y);
2122     FREE(t);
2123         sciprint(_("%s: No more memory.\n"),"multifrontal_supernodal_solve_llt");
2124     return -1;
2125   }
2126
2127   for (i=0; i<L->n; i++) x[i] = b[i];
2128
2129   recursive_supernodal_solve_l (L->n_sn,
2130                                 TRUE,  /* this is the root */
2131                                 L->first_child, L->next_child,
2132                                 L->sn_struct,L->sn_size,L->sn_up_size,
2133                                 L->sn_blocks_ld, L->sn_blocks,
2134                                 L->up_blocks_ld, L->up_blocks,
2135                                 y, x, t);
2136
2137   recursive_supernodal_solve_lt(L->n_sn,
2138                                 TRUE,  /* this is the root */
2139                                 L->first_child, L->next_child,
2140                                 L->sn_struct,L->sn_size,L->sn_up_size,
2141                                 L->sn_blocks_ld, L->sn_blocks,
2142                                 L->up_blocks_ld, L->up_blocks,
2143                                 x, y, t);
2144
2145   FREE(y);
2146   FREE(t);
2147     
2148   return 0;
2149 }
2150
2151 /*****************************************************
2152  *  some utility routines                            *
2153  *****************************************************/
2154
2155 taucs_ccs_matrix*
2156 taucs_supernodal_factor_to_ccs(void* vL)
2157 {
2158   supernodal_factor_matrix* L = (supernodal_factor_matrix*) vL;
2159   taucs_ccs_matrix* C;
2160   int n,nnz;
2161   int i,j,ip,jp,sn,next;
2162   double v;
2163   int* len;
2164
2165   n = L->n;
2166
2167   len = (int*) MALLOC(n*sizeof(int));
2168   if (!len) return NULL;
2169
2170   nnz = 0;
2171
2172   for (sn=0; sn<L->n_sn; sn++) {
2173     for (jp=0; jp<(L->sn_size)[sn]; jp++) {
2174       j = (L->sn_struct)[sn][jp];
2175       len[j] = 0;
2176
2177       for (ip=jp; ip<(L->sn_size)[sn]; ip++) {
2178         v = (L->sn_blocks)[sn][ jp*(L->sn_blocks_ld)[sn] + ip ];
2179
2180         if (v) { 
2181           len[j] ++;
2182           nnz ++;
2183         }
2184       }
2185       for (ip=(L->sn_size)[sn]; ip<(L->sn_up_size)[sn]; ip++) {
2186         v = (L->up_blocks)[sn][ jp*(L->up_blocks_ld)[sn] + (ip-(L->sn_size)[sn]) ];
2187
2188         if (v) { 
2189           len[j] ++;
2190           nnz ++;
2191         }
2192       }
2193     }
2194   }
2195
2196
2197   C = taucs_ccs_create(n,n,nnz);
2198   if (!C) {
2199     FREE(len);
2200     return NULL;
2201   }
2202   C->flags = TAUCS_TRIANGULAR | TAUCS_LOWER;
2203
2204   (C->colptr)[0] = 0;
2205   for (j=1; j<=n; j++) (C->colptr)[j] = (C->colptr)[j-1] + len[j-1];
2206
2207   FREE(len);
2208
2209   for (sn=0; sn<L->n_sn; sn++) {
2210     for (jp=0; jp<(L->sn_size)[sn]; jp++) {
2211       j = (L->sn_struct)[sn][jp];
2212
2213       next = (C->colptr)[j];
2214
2215       for (ip=jp; ip<(L->sn_size)[sn]; ip++) {
2216         i = (L->sn_struct)[sn][ ip ];
2217         v = (L->sn_blocks)[sn][ jp*(L->sn_blocks_ld)[sn] + ip ];
2218
2219         if (v == 0.0) continue;
2220
2221         (C->rowind)[next] = i;
2222         (C->values)[next] = v;
2223         next++;
2224       }
2225       for (ip=(L->sn_size)[sn]; ip<(L->sn_up_size)[sn]; ip++) {
2226         i = (L->sn_struct)[sn][ ip ];
2227         v = (L->up_blocks)[sn][ jp*(L->up_blocks_ld)[sn] + (ip-(L->sn_size)[sn]) ];
2228
2229         if (v == 0.0) continue;
2230
2231         (C->rowind)[next] = i;
2232         (C->values)[next] = v;
2233         next++;
2234       }
2235     }
2236   }
2237
2238   return C;
2239 }
2240
2241 /* 
2242  *   from the preeeding routine here is a small one to get only
2243  *   the number of non zero elements of the factor (Bruno le 29/11/2001)
2244  *   In fact the supernodal struct keeps on some exact zeros and so the
2245  *   nnz get here is representative of the memory used by the supernodal
2246  *   struc.
2247  */
2248
2249 int taucs_get_nnz_from_supernodal_factor(void* vL)
2250 {
2251   supernodal_factor_matrix* L = (supernodal_factor_matrix*) vL;
2252   int nnz = 0;
2253   int jp = 0, sn = 0;
2254
2255   nnz = 0;
2256
2257   for (sn=0; sn<L->n_sn; sn++)
2258     for (jp=0; jp<(L->sn_size)[sn]; jp++) 
2259       nnz += (L->sn_up_size)[sn] - jp;
2260
2261   return ( nnz );
2262 }
2263
2264
2265