dummy/solver/ascgauss.c

/*
 *  ascgauss: sparse gauss elimination
 *  by Ken Tyner and Ben Allan
 *  Created: 4/14/97
 *  Version: $Revision: 1.6 $
 *  Version control file: $RCSfile: ascgauss.c,v $
 *  Date last modified: $Date: 2000/01/25 02:26:48 $
 *  Last modified by: $Author: ballan $
 *
 *  This file is part of the SLV solver.
 *  The SLV solver is free software; you can redistribute
 *  it and/or modify it under the terms of the GNU General Public License as
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version.
 *
 *  The SLV solver is distributed in hope that it will be
 *  useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along with
 *  the program; if not, write to the Free Software Foundation, Inc., 675
 *  Mass Ave, Cambridge, MA 02139 USA.  Check the file named COPYING.
 *  COPYING is found in ../compiler.
 */

/*
 * This file contains gauss and modified gauss implementations.
 * The first gauss is a partial column pivoted gauss L\U.
 * Initially (4/97) it is not the fastest possible code, but is
 * to be easily understood.
 * It needs a header.
 */
#define GAUSSALONE 0
#if GAUSSALONE
#include <math.h>
#include <utilities/ascConfig.h>
#include <compiler/compiler.h>
#include <utilities/ascMalloc.h>
#include <utilities/mem.h>
#include <general/tm_time.h>
#include "mtx.h"
#include "linsolqr.h"
#define D_ZERO (real64)0.0
#define ZERO (int32)0
#define D_ONE (real64)1.0
#endif

#define GAUSS_DEBUG TRUE

/*
 * See the comments attached to forward_substitute.
 * This code otherwise is the same. It uses a 2 bodied
 * matrix as well as making use of mtx_ALL_COLS and
 * mtx_ALL_ROWS whereever possible.
 * We make the following additional assumptions here:
 *   (if they do not hold, do NOT use this function)
 * - sys->inverse (L) has no incidence that is not multipliers
 *   (and the diagonal of 1s is NOT in sys->inverse.)
 * - sys->factors (U) has no incidence on the upper triangle,
 *   including the diagonal, or outside the factored region.
 *    relaxation: incidence anywhere allowed if value = 0.0
 *                since 0 doesn't contribute to a dot product
 *                and the only thing we do with triangles is dot them.
 * - There may be singular rows and columns in the factorization,
 *   but any additions coming from these rows/columns during
 *   mtx_ALL_*O*S operations will not contribute to sums because the
 *   user zeroed the arr entries corresponding to these before
 *   calling this function.
 */
static void forward_substitute_gauss2(linsolqr_system_t sys,
                     real64 *arr,
                     boolean transpose)
{
  mtx_coord_t nz;
  real64 sum, *pivlist;
  mtx_matrix_t mtx;
  int32 dotlim;
  boolean nonzero_found=FALSE;

  pivlist=sys->ludata->pivlist;
  dotlim = sys->rng.low+sys->rank;
  if (transpose) {      /* arr is indexed by original column number */
    mtx=sys->factors;
    for( nz.col=sys->rng.low; nz.col < dotlim; ++(nz.col) ) {
      register int32 org_col;

      org_col = mtx_col_to_org(mtx,nz.col);
      if (arr[org_col]!=D_ZERO) nonzero_found=TRUE;
      if (nonzero_found) {
    sum=mtx_col_dot_full_org_vec(mtx,nz.col,arr,mtx_ALL_ROWS,TRUE);
    arr[org_col] = (arr[org_col] - sum) / pivlist[nz.col];
      }
    }
  } else {          /* arr is indexed by original row number */
    mtx=sys->inverse;
    for( nz.row=sys->rng.low; nz.row < dotlim; ++(nz.row) ) {
      register int32 org_row;

      org_row = mtx_row_to_org(mtx,nz.row);
      if (arr[org_row]!=D_ZERO) nonzero_found=TRUE;
      if (nonzero_found) {
    sum = mtx_row_dot_full_org_vec(mtx,nz.row,arr,mtx_ALL_COLS,TRUE);
    arr[org_row] -= sum;
      }
    }
  }
}

/*
 * See the comments attached to backward_substitute and
 * forward_substitute2.
 * When solving for the transpose, then we are actually
 * running of the lower triangle, hence we use sys->factors.
 * Otherwise we use sys->inverse which stores U.
 */
static void backward_substitute_gauss2(linsolqr_system_t sys,
                      real64 *arr,
                      boolean transpose)
{
  mtx_coord_t nz;
  real64 sum, *pivlist;
  mtx_matrix_t mtx;
  int32 dotlim;
  boolean nonzero_found=FALSE;  /* once TRUE, substitution must be done
                   over remaining rows/cols */

  dotlim=sys->rng.low;
  pivlist=sys->ludata->pivlist;
  if (transpose) {      /* arr is indexed by original column number */
    mtx = sys->inverse;
    for( nz.col = sys->rng.low+sys->rank-1; nz.col >= dotlim ; --(nz.col) ) {
      register int32 org_col;

      org_col = mtx_col_to_org(mtx,nz.col);
      if (arr[org_col] != D_ZERO) nonzero_found=TRUE;
      if (nonzero_found) {
    sum = mtx_col_dot_full_org_vec(mtx,nz.col,arr,mtx_ALL_ROWS,TRUE);
    arr[org_col] -= sum;
      }
    }
  } else {          /* arr is indexed by original row number */
    /* we are working from the bottom up */
    mtx = sys->factors;
    for( nz.row = sys->rng.low+sys->rank-1; nz.row >= dotlim ; --(nz.row) ) {
      register int32 org_row;

      org_row = mtx_row_to_org(mtx,nz.row);
      if (arr[org_row]!=D_ZERO) nonzero_found=TRUE;
      if (nonzero_found) {
    sum= mtx_row_dot_full_org_vec(mtx,nz.row,arr,mtx_ALL_COLS,TRUE);
    arr[org_row] = (arr[org_row] - sum) / pivlist[nz.row];
      }
    }
  }
}

static int gauss2_solve(linsolqr_system_t sys, struct rhs_list *rl)
{
  /* zero any unsolved for vars first so they don't contaminate
     mtx_ALL_*O*S dot products.
  */
  zero_unpivoted_vars(sys,rl->varvalue,rl->transpose);
  forward_substitute_gauss2(sys,rl->varvalue,rl->transpose);
  backward_substitute_gauss2(sys,rl->varvalue,rl->transpose);
  return 0;
}


/* this function eliminates a column */
static
int eliminate_column2(
mtx_matrix_t mtx, /* the matrix being factored (upper triangularized) */
mtx_matrix_t lmtx,  /* the multipliers of the factorization, L */
int32 column,   /* the column being eliminated (the pivot is (column,column))*/
int32 lastrow,  /* rows from column+1 to lastrow will be eliminated in column*/
real64 pivot,   /* the pivot value */
mtx_sparse_t *pivotrow  /* the pivot row, less the pivot entry, orgindexed */
)
{
  mtx_coord_t nz;
  mtx_range_t elimrng;
  real64 multiplier;

  /* for rows in column+1 .. lastrow */
  elimrng.low = column+1;
  elimrng.high = lastrow;

  nz.row = mtx_FIRST;
  while ((void)mtx_next_in_col(mtx,&nz,&elimrng),
         nz.row !=mtx_LAST) {
    multiplier = mtx_value(mtx,&nz)/pivot;
    mtx_add_row_sparse(mtx,nz.row,-multiplier,pivotrow);
    mtx_fill_value(lmtx,&nz,multiplier);
  }
  mtx_clear_col(mtx,column,&elimrng);
  return 0;
}

/* This factors a square region given in sys->rng.
 * We might like to have a wider pivot region later, so does not assume
 * that the region to the right of the factorization is empty.
 * Assumes the region to the left of the factorization IS empty.
 */
int gaussba2_factor_square(linsolqr_system_t sys)
{
  mtx_coord_t nz;
  int32 lastrow,lastcolumn,column;
  mtx_range_t pivot_candidates;
  real64 *tmp;
  real64 pivot, *pivots;
  int32 length, defect = 0;
  mtx_matrix_t mtx, lmtx;
  mtx_sparse_t *pivotrow;
  int err = 0;

  length = sys->rng.high - sys->rng.low + 1; /* ? */
  sys->smallest_pivot = MAXDOUBLE;
  lastrow = lastcolumn = pivot_candidates.high = sys->rng.high;
  pivot_candidates.low = sys->rng.low;

  tmp = sys->ludata->tmp;
  pivots = sys->ludata->pivlist;
  mtx = sys->factors;
  lmtx = sys->inverse;
  pivotrow = sys->ludata->tmp_sparse;

/* need pivotrow allocation */

  for (column = sys->rng.low; column <= lastcolumn; ) {
    /* pick pivot. */
    nz.row = column;
    pivot = mtx_get_pivot_col(mtx,&nz,&pivot_candidates,&(pivots[nz.row]),
                              sys->ptol,sys->pivot_zero);

    if (  pivot < sys->pivot_zero ) { /* pivot is an epsilon */
      /* If no pivot, drag to bottom and lastcolumn-- */
      mtx_drag(mtx, nz.row, lastcolumn);
      number_drag(pivots, nz.row, lastcolumn);
      lastcolumn--;
      pivot_candidates.high--; /* is this correct? */
      defect++;
    } else {
      /* collect pivot row. pivot removed to reduce numerics cost*/
      mtx_clear_coord(mtx,nz.row,nz.col);
      pivotrow = mtx_org_row_sparse(mtx,nz.row,pivotrow,
                    mtx_ALL_COLS,mtx_IGNORE_ZEROES); /* check this */
      /* eliminate column */
#if (GAUSS_DEBUG)
  {
    FILE *fp;
    char fname[80];

    sprintf(fname,"/tmp/lu/uga.%d",column");
    fp = fopen(fname,"w+");
    mtx_write_region_human_rows(fp,mtx,mtx_ENTIRE_MATRIX);
    fclose(fp);

    sprintf(fname,"/tmp/lu/lga.%d",column");
    fp = fopen(fname,"w+");
    mtx_write_region_human_rows(fp,lmtx,mtx_ENTIRE_MATRIX);
    fclose(fp);

    sprintf(fname,"/tmp/lu/aga.plot%d",column);
    fp=fopen(fname,"w+");
    mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
    fclose(fp);
  }
#endif
      err = eliminate_column2(mtx,lmtx,column,lastrow,pivot,pivotrow);
      if (err != 0) {
        break;
      }
      pivot_candidates.low++;
      column++;
    }
  }
  sys->rank = pivot_candidates.high - sys->rng.low + 1;
  return err; /* err = 0 is good */
}


int gauss2_entry(linsolqr_system_t sys, mtx_region_t *region)
{
  struct rhs_list *rl;
  double time;

  CHECK_SYSTEM(sys);
  if (sys->factored) return 0;
  switch(sys->fmethod) {
  case gauss_ba2:    /* add new methods here */
    break;
  default:
    return 1;
  }

  if (ISNULL(sys->ludata)) return 1;
  if (NOTNULL(sys->inverse)) mtx_destroy(sys->inverse);
  sys->inverse = NULL;
  if (NOTNULL(sys->factors)) mtx_destroy(sys->factors);
  if (region == mtx_ENTIRE_MATRIX) determine_pivot_range(sys);
  else square_region(sys,region);

  sys->factors = mtx_copy_region(sys->coef,region);
  sys->inverse = mtx_create_slave(sys->factors);
  sys->rank = -1;
  sys->smallest_pivot = MAXDOUBLE;
  for (rl = sys->rl ; NOTNULL(rl) ; rl = rl->next)
    rl->solved = FALSE;
  insure_capacity(sys);
  insure_lu_capacity(sys);

  time = tm_cpu_time();
  switch(sys->fmethod) {
  case gauss_ba2:
    gaussba2_factor_square(sys);
    break;
  default:
    return 1;
  }
  sys->factored = TRUE;

#define KAA_DEBUG
#ifdef KAA_DEBUG
  time = tm_cpu_time() - time;
  FPRINTF(stderr,"Time for Inversion = %f\n",time);
  FPRINTF(stderr,"Non-zeros in A = %d\n",
    mtx_nonzeros_in_region(sys->coef,region));
  FPRINTF(stderr,"Non-zeros in U+L = %d\n",
    mtx_nonzeros_in_region(sys->factors,region) +
    mtx_nonzeros_in_region(sys->inverse,0));
#endif /* KAA_DEBUG */
#undef KAA_DEBUG
  return 0;
}


/* The following crap is from kirk. we should be ditching it and writing a
 * real gauss suitable to ken's and other needs.
 */
#define GAUSS_ELIMINATE_ROWS_IMPLEMENTED 0
#if GAUSS_ELIMINATE_ROWS_IMPLEMENTED
static void gauss_eliminate_rows(mtx_matrix_t mtx,
                 mtx_matrix_t lower_mtx,
                 int32 col,
                 real64 *value,
                 int32 *index,
                 real64 pivot)
{
  mtx_coord_t nz;
  real64 multiplier;
  int32 i, ii;
  int colcount, row;

/*
  colcount = mtx_cur_col_to_csr(mtx,col,value,index);
*/
  mtx_clear_col(mtx,col,mtx_ALL_ROWS);

  row = col;
  nz.col = col;
  for (ii=0; ii<colcount; ii++) {
    i = index[ii];
    if (i <= row) {
      /* nothing to do - we would be back patching U here. */
      continue;
    }
    else{
      multiplier = value[ii]/pivot;
      mtx_add_row(mtx,row,i,-multiplier,mtx_ALL_COLS);
      nz.row = i;
      mtx_fill_value(lower_mtx,&nz,multiplier);
    }
  }
}


/*
 * At the moment this code assumes that we are given a
 * square region. It also assumes that the submatrix A11 will
 * be provided; in which case pivot selection should be restricted
 * to A11. Later on we will modify the code to accept a rectangular
 * region, we will then make use of the 'last_pivotable column',
 * and sys->reg. At the moment only sys->rng and A11 are used.
 */
static int kirk_gauss_factor(linsolqr_system_t sys,mtx_region_t *A11)
{
  /*
   * THIS HAS BUG -- row can not be updated.
   */
  mtx_matrix_t mtx, lower_mtx;
  mtx_coord_t nz;
  int32 last_col;
  mtx_range_t pivot_candidates;
  real64 pivot, abs_pivot, *pivots, max;
  int32 row,col;
  int32 *index = NULL;      /* used by csr */
  real64 *value = NULL;     /* used by csr */
  int length;

  mtx = sys->factors;
  lower_mtx = sys->inverse;
  pivots = sys->ludata->pivlist;

  last_col = sys->rng.high;
  pivot_candidates.high = A11->col.high;    /* initialize */

  length = sys->rng.high - sys->rng.low + 1;
  value = (real64 *)ascmalloc(length*sizeof(real64));
  index = (int32 *)ascmalloc(length*sizeof(int32));

  /*
   * First find the maximum element in the col *regardless*. This
   * will leave nz at the correct position in case we need to a
   * column swap. We have to adjust the pivot range as we go along.
   */
  for (col = sys->rng.low; col < last_col; ) {
    nz.row = nz.col = col;
    pivot_candidates.low = col;
    if (nz.col > A11->col.high) {
      pivot_candidates.high = sys->rng.high;
    }
    pivot = mtx_value(mtx,&nz);

    max = mtx_col_max(mtx, &nz, &pivot_candidates, &pivots[col]);
    if (max <= sys->pivot_zero) {
      FPRINTF(stderr,"Matrix appears numerically singular at %d %g %g\n",
          row,pivot,max);
      /* patch for the time being */
      pivot = pivots[col] = 1.0;
      col++;
      continue;
    }

    abs_pivot = fabs(pivot);
    if ((abs_pivot >= sys->ptol * max)) {
      pivots[col] = pivot;
    }
    else {
      mtx_swap_rows(mtx, row, nz.row);
      pivot = pivots[col] = max;
    }

    gauss_eliminate_rows(mtx, lower_mtx, col, value, index, pivot);
    col++;
  }

  if (value) ascfree((char *)value);
  if (index) ascfree((char *)index);

  sys->rank = last_col - sys->rng.low + 1;
  return sys->rank;
}
#endif /* GAUSS_ELIMINATE_ROWS_IMPLEMENTED */
1	/*
2	* ascgauss: sparse gauss elimination
3	* by Ken Tyner and Ben Allan
4	* Created: 4/14/97
5	* Version: $Revision: 1.6 $
6	* Version control file: $RCSfile: ascgauss.c,v $
7	* Date last modified: $Date: 2000/01/25 02:26:48 $
8	* Last modified by: $Author: ballan $
9	*
10	* This file is part of the SLV solver.
11	* The SLV solver is free software; you can redistribute
12	* it and/or modify it under the terms of the GNU General Public License as
13	* published by the Free Software Foundation; either version 2 of the
14	* License, or (at your option) any later version.
15	*
16	* The SLV solver is distributed in hope that it will be
17	* useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
18	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19	* General Public License for more details.
20	*
21	* You should have received a copy of the GNU General Public License along with
22	* the program; if not, write to the Free Software Foundation, Inc., 675
23	* Mass Ave, Cambridge, MA 02139 USA. Check the file named COPYING.
24	* COPYING is found in ../compiler.
25	*/
26
27	/*
28	* This file contains gauss and modified gauss implementations.
29	* The first gauss is a partial column pivoted gauss L\U.
30	* Initially (4/97) it is not the fastest possible code, but is
31	* to be easily understood.
32	* It needs a header.
33	*/
34	#define GAUSSALONE 0
35	#if GAUSSALONE
36	#include <math.h>
37	#include <utilities/ascConfig.h>
38	#include <compiler/compiler.h>
39	#include <utilities/ascMalloc.h>
40	#include <utilities/mem.h>
41	#include <general/tm_time.h>
42	#include "mtx.h"
43	#include "linsolqr.h"
44	#define D_ZERO (real64)0.0
45	#define ZERO (int32)0
46	#define D_ONE (real64)1.0
47	#endif
48
49	#define GAUSS_DEBUG TRUE
50
51	/*
52	* See the comments attached to forward_substitute.
53	* This code otherwise is the same. It uses a 2 bodied
54	* matrix as well as making use of mtx_ALL_COLS and
55	* mtx_ALL_ROWS whereever possible.
56	* We make the following additional assumptions here:
57	* (if they do not hold, do NOT use this function)
58	* - sys->inverse (L) has no incidence that is not multipliers
59	* (and the diagonal of 1s is NOT in sys->inverse.)
60	* - sys->factors (U) has no incidence on the upper triangle,
61	* including the diagonal, or outside the factored region.
62	* relaxation: incidence anywhere allowed if value = 0.0
63	* since 0 doesn't contribute to a dot product
64	* and the only thing we do with triangles is dot them.
65	* - There may be singular rows and columns in the factorization,
66	* but any additions coming from these rows/columns during
67	* mtx_ALL_OS operations will not contribute to sums because the
68	* user zeroed the arr entries corresponding to these before
69	* calling this function.
70	*/
71	static void forward_substitute_gauss2(linsolqr_system_t sys,
72	real64 *arr,
73	boolean transpose)
74	{
75	mtx_coord_t nz;
76	real64 sum, *pivlist;
77	mtx_matrix_t mtx;
78	int32 dotlim;
79	boolean nonzero_found=FALSE;
80
81	pivlist=sys->ludata->pivlist;
82	dotlim = sys->rng.low+sys->rank;
83	if (transpose) { /* arr is indexed by original column number */
84	mtx=sys->factors;
85	for( nz.col=sys->rng.low; nz.col < dotlim; ++(nz.col) ) {
86	register int32 org_col;
87
88	org_col = mtx_col_to_org(mtx,nz.col);
89	if (arr[org_col]!=D_ZERO) nonzero_found=TRUE;
90	if (nonzero_found) {
91	sum=mtx_col_dot_full_org_vec(mtx,nz.col,arr,mtx_ALL_ROWS,TRUE);
92	arr[org_col] = (arr[org_col] - sum) / pivlist[nz.col];
93	}
94	}
95	} else { /* arr is indexed by original row number */
96	mtx=sys->inverse;
97	for( nz.row=sys->rng.low; nz.row < dotlim; ++(nz.row) ) {
98	register int32 org_row;
99
100	org_row = mtx_row_to_org(mtx,nz.row);
101	if (arr[org_row]!=D_ZERO) nonzero_found=TRUE;
102	if (nonzero_found) {
103	sum = mtx_row_dot_full_org_vec(mtx,nz.row,arr,mtx_ALL_COLS,TRUE);
104	arr[org_row] -= sum;
105	}
106	}
107	}
108	}
109
110	/*
111	* See the comments attached to backward_substitute and
112	* forward_substitute2.
113	* When solving for the transpose, then we are actually
114	* running of the lower triangle, hence we use sys->factors.
115	* Otherwise we use sys->inverse which stores U.
116	*/
117	static void backward_substitute_gauss2(linsolqr_system_t sys,
118	real64 *arr,
119	boolean transpose)
120	{
121	mtx_coord_t nz;
122	real64 sum, *pivlist;
123	mtx_matrix_t mtx;
124	int32 dotlim;
125	boolean nonzero_found=FALSE; /* once TRUE, substitution must be done
126	over remaining rows/cols */
127
128	dotlim=sys->rng.low;
129	pivlist=sys->ludata->pivlist;
130	if (transpose) { /* arr is indexed by original column number */
131	mtx = sys->inverse;
132	for( nz.col = sys->rng.low+sys->rank-1; nz.col >= dotlim ; --(nz.col) ) {
133	register int32 org_col;
134
135	org_col = mtx_col_to_org(mtx,nz.col);
136	if (arr[org_col] != D_ZERO) nonzero_found=TRUE;
137	if (nonzero_found) {
138	sum = mtx_col_dot_full_org_vec(mtx,nz.col,arr,mtx_ALL_ROWS,TRUE);
139	arr[org_col] -= sum;
140	}
141	}
142	} else { /* arr is indexed by original row number */
143	/* we are working from the bottom up */
144	mtx = sys->factors;
145	for( nz.row = sys->rng.low+sys->rank-1; nz.row >= dotlim ; --(nz.row) ) {
146	register int32 org_row;
147
148	org_row = mtx_row_to_org(mtx,nz.row);
149	if (arr[org_row]!=D_ZERO) nonzero_found=TRUE;
150	if (nonzero_found) {
151	sum= mtx_row_dot_full_org_vec(mtx,nz.row,arr,mtx_ALL_COLS,TRUE);
152	arr[org_row] = (arr[org_row] - sum) / pivlist[nz.row];
153	}
154	}
155	}
156	}
157
158	static int gauss2_solve(linsolqr_system_t sys, struct rhs_list *rl)
159	{
160	/* zero any unsolved for vars first so they don't contaminate
161	mtx_ALL_OS dot products.
162	*/
163	zero_unpivoted_vars(sys,rl->varvalue,rl->transpose);
164	forward_substitute_gauss2(sys,rl->varvalue,rl->transpose);
165	backward_substitute_gauss2(sys,rl->varvalue,rl->transpose);
166	return 0;
167	}
168
169
170
171	/* this function eliminates a column */
172	static
173	int eliminate_column2(
174	mtx_matrix_t mtx, /* the matrix being factored (upper triangularized) */
175	mtx_matrix_t lmtx, /* the multipliers of the factorization, L */
176	int32 column, /* the column being eliminated (the pivot is (column,column))*/
177	int32 lastrow, /* rows from column+1 to lastrow will be eliminated in column*/
178	real64 pivot, /* the pivot value */
179	mtx_sparse_t pivotrow / the pivot row, less the pivot entry, orgindexed */
180	)
181	{
182	mtx_coord_t nz;
183	mtx_range_t elimrng;
184	real64 multiplier;
185
186	/* for rows in column+1 .. lastrow */
187	elimrng.low = column+1;
188	elimrng.high = lastrow;
189
190	nz.row = mtx_FIRST;
191	while ((void)mtx_next_in_col(mtx,&nz,&elimrng),
192	nz.row !=mtx_LAST) {
193	multiplier = mtx_value(mtx,&nz)/pivot;
194	mtx_add_row_sparse(mtx,nz.row,-multiplier,pivotrow);
195	mtx_fill_value(lmtx,&nz,multiplier);
196	}
197	mtx_clear_col(mtx,column,&elimrng);
198	return 0;
199	}
200
201	/* This factors a square region given in sys->rng.
202	* We might like to have a wider pivot region later, so does not assume
203	* that the region to the right of the factorization is empty.
204	* Assumes the region to the left of the factorization IS empty.
205	*/
206	int gaussba2_factor_square(linsolqr_system_t sys)
207	{
208	mtx_coord_t nz;
209	int32 lastrow,lastcolumn,column;
210	mtx_range_t pivot_candidates;
211	real64 *tmp;
212	real64 pivot, *pivots;
213	int32 length, defect = 0;
214	mtx_matrix_t mtx, lmtx;
215	mtx_sparse_t *pivotrow;
216	int err = 0;
217
218	length = sys->rng.high - sys->rng.low + 1; /* ? */
219	sys->smallest_pivot = MAXDOUBLE;
220	lastrow = lastcolumn = pivot_candidates.high = sys->rng.high;
221	pivot_candidates.low = sys->rng.low;
222
223	tmp = sys->ludata->tmp;
224	pivots = sys->ludata->pivlist;
225	mtx = sys->factors;
226	lmtx = sys->inverse;
227	pivotrow = sys->ludata->tmp_sparse;
228
229	/* need pivotrow allocation */
230
231	for (column = sys->rng.low; column <= lastcolumn; ) {
232	/* pick pivot. */
233	nz.row = column;
234	pivot = mtx_get_pivot_col(mtx,&nz,&pivot_candidates,&(pivots[nz.row]),
235	sys->ptol,sys->pivot_zero);
236
237	if ( pivot < sys->pivot_zero ) { /* pivot is an epsilon */
238	/* If no pivot, drag to bottom and lastcolumn-- */
239	mtx_drag(mtx, nz.row, lastcolumn);
240	number_drag(pivots, nz.row, lastcolumn);
241	lastcolumn--;
242	pivot_candidates.high--; /* is this correct? */
243	defect++;
244	} else {
245	/* collect pivot row. pivot removed to reduce numerics cost*/
246	mtx_clear_coord(mtx,nz.row,nz.col);
247	pivotrow = mtx_org_row_sparse(mtx,nz.row,pivotrow,
248	mtx_ALL_COLS,mtx_IGNORE_ZEROES); /* check this */
249	/* eliminate column */
250	#if (GAUSS_DEBUG)
251	{
252	FILE *fp;
253	char fname[80];
254
255	sprintf(fname,"/tmp/lu/uga.%d",column");
256	fp = fopen(fname,"w+");
257	mtx_write_region_human_rows(fp,mtx,mtx_ENTIRE_MATRIX);
258	fclose(fp);
259
260	sprintf(fname,"/tmp/lu/lga.%d",column");
261	fp = fopen(fname,"w+");
262	mtx_write_region_human_rows(fp,lmtx,mtx_ENTIRE_MATRIX);
263	fclose(fp);
264
265	sprintf(fname,"/tmp/lu/aga.plot%d",column);
266	fp=fopen(fname,"w+");
267	mtx_write_region_plot(fp,mtx,mtx_ENTIRE_MATRIX);
268	fclose(fp);
269	}
270	#endif
271	err = eliminate_column2(mtx,lmtx,column,lastrow,pivot,pivotrow);
272	if (err != 0) {
273	break;
274	}
275	pivot_candidates.low++;
276	column++;
277	}
278	}
279	sys->rank = pivot_candidates.high - sys->rng.low + 1;
280	return err; /* err = 0 is good */
281	}
282
283
284	int gauss2_entry(linsolqr_system_t sys, mtx_region_t *region)
285	{
286	struct rhs_list *rl;
287	double time;
288
289	CHECK_SYSTEM(sys);
290	if (sys->factored) return 0;
291	switch(sys->fmethod) {
292	case gauss_ba2: /* add new methods here */
293	break;
294	default:
295	return 1;
296	}
297
298	if (ISNULL(sys->ludata)) return 1;
299	if (NOTNULL(sys->inverse)) mtx_destroy(sys->inverse);
300	sys->inverse = NULL;
301	if (NOTNULL(sys->factors)) mtx_destroy(sys->factors);
302	if (region == mtx_ENTIRE_MATRIX) determine_pivot_range(sys);
303	else square_region(sys,region);
304
305	sys->factors = mtx_copy_region(sys->coef,region);
306	sys->inverse = mtx_create_slave(sys->factors);
307	sys->rank = -1;
308	sys->smallest_pivot = MAXDOUBLE;
309	for (rl = sys->rl ; NOTNULL(rl) ; rl = rl->next)
310	rl->solved = FALSE;
311	insure_capacity(sys);
312	insure_lu_capacity(sys);
313
314	time = tm_cpu_time();
315	switch(sys->fmethod) {
316	case gauss_ba2:
317	gaussba2_factor_square(sys);
318	break;
319	default:
320	return 1;
321	}
322	sys->factored = TRUE;
323
324	#define KAA_DEBUG
325	#ifdef KAA_DEBUG
326	time = tm_cpu_time() - time;
327	FPRINTF(stderr,"Time for Inversion = %f\n",time);
328	FPRINTF(stderr,"Non-zeros in A = %d\n",
329	mtx_nonzeros_in_region(sys->coef,region));
330	FPRINTF(stderr,"Non-zeros in U+L = %d\n",
331	mtx_nonzeros_in_region(sys->factors,region) +
332	mtx_nonzeros_in_region(sys->inverse,0));
333	#endif /* KAA_DEBUG */
334	#undef KAA_DEBUG
335	return 0;
336	}
337
338
339
340
341
342
343
344
345
346	/* The following crap is from kirk. we should be ditching it and writing a
347	* real gauss suitable to ken's and other needs.
348	*/
349	#define GAUSS_ELIMINATE_ROWS_IMPLEMENTED 0
350	#if GAUSS_ELIMINATE_ROWS_IMPLEMENTED
351	static void gauss_eliminate_rows(mtx_matrix_t mtx,
352	mtx_matrix_t lower_mtx,
353	int32 col,
354	real64 *value,
355	int32 *index,
356	real64 pivot)
357	{
358	mtx_coord_t nz;
359	real64 multiplier;
360	int32 i, ii;
361	int colcount, row;
362
363	/*
364	colcount = mtx_cur_col_to_csr(mtx,col,value,index);
365	*/
366	mtx_clear_col(mtx,col,mtx_ALL_ROWS);
367
368	row = col;
369	nz.col = col;
370	for (ii=0; ii<colcount; ii++) {
371	i = index[ii];
372	if (i <= row) {
373	/* nothing to do - we would be back patching U here. */
374	continue;
375	}
376	else{
377	multiplier = value[ii]/pivot;
378	mtx_add_row(mtx,row,i,-multiplier,mtx_ALL_COLS);
379	nz.row = i;
380	mtx_fill_value(lower_mtx,&nz,multiplier);
381	}
382	}
383	}
384
385
386	/*
387	* At the moment this code assumes that we are given a
388	* square region. It also assumes that the submatrix A11 will
389	* be provided; in which case pivot selection should be restricted
390	* to A11. Later on we will modify the code to accept a rectangular
391	* region, we will then make use of the 'last_pivotable column',
392	* and sys->reg. At the moment only sys->rng and A11 are used.
393	*/
394	static int kirk_gauss_factor(linsolqr_system_t sys,mtx_region_t *A11)
395	{
396	/*
397	* THIS HAS BUG -- row can not be updated.
398	*/
399	mtx_matrix_t mtx, lower_mtx;
400	mtx_coord_t nz;
401	int32 last_col;
402	mtx_range_t pivot_candidates;
403	real64 pivot, abs_pivot, *pivots, max;
404	int32 row,col;
405	int32 index = NULL; / used by csr */
406	real64 value = NULL; / used by csr */
407	int length;
408
409	mtx = sys->factors;
410	lower_mtx = sys->inverse;
411	pivots = sys->ludata->pivlist;
412
413	last_col = sys->rng.high;
414	pivot_candidates.high = A11->col.high; /* initialize */
415
416	length = sys->rng.high - sys->rng.low + 1;
417	value = (real64 )ascmalloc(lengthsizeof(real64));
418	index = (int32 )ascmalloc(lengthsizeof(int32));
419
420	/*
421	* First find the maximum element in the col regardless. This
422	* will leave nz at the correct position in case we need to a
423	* column swap. We have to adjust the pivot range as we go along.
424	*/
425	for (col = sys->rng.low; col < last_col; ) {
426	nz.row = nz.col = col;
427	pivot_candidates.low = col;
428	if (nz.col > A11->col.high) {
429	pivot_candidates.high = sys->rng.high;
430	}
431	pivot = mtx_value(mtx,&nz);
432
433	max = mtx_col_max(mtx, &nz, &pivot_candidates, &pivots[col]);
434	if (max <= sys->pivot_zero) {
435	FPRINTF(stderr,"Matrix appears numerically singular at %d %g %g\n",
436	row,pivot,max);
437	/* patch for the time being */
438	pivot = pivots[col] = 1.0;
439	col++;
440	continue;
441	}
442
443	abs_pivot = fabs(pivot);
444	if ((abs_pivot >= sys->ptol * max)) {
445	pivots[col] = pivot;
446	}
447	else {
448	mtx_swap_rows(mtx, row, nz.row);
449	pivot = pivots[col] = max;
450	}
451
452	gauss_eliminate_rows(mtx, lower_mtx, col, value, index, pivot);
453	col++;
454	}
455
456	if (value) ascfree((char *)value);
457	if (index) ascfree((char *)index);
458
459	sys->rank = last_col - sys->rng.low + 1;
460	return sys->rank;
461	}
462	#endif /* GAUSS_ELIMINATE_ROWS_IMPLEMENTED */