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 "solver/mtx.h"
#include "solver/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	aw0a	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 "solver/mtx.h"
43			#include "solver/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 */