generic/interface/Sensitivity.c

/*
 *  Sensitivity.c
 *  by Kirk Abbott and Ben Allan
 *  Created: 1/94
 *  Version: $Revision: 1.31 $
 *  Version control file: $RCSfile: Sensitivity.c,v $
 *  Date last modified: $Date: 2003/08/23 18:43:08 $
 *  Last modified by: $Author: ballan $
 *
 *  This file is part of the ASCEND Tcl/Tk interface
 *
 *  Copyright 1997, Carnegie Mellon University
 *
 *  The ASCEND Tcl/Tk interface 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 ASCEND Tcl/Tk interface 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.
 */


/*
 *  Sensititvity analysis code.
 */

#include <math.h>
#include <tcl.h>
#include <utilities/ascConfig.h>
#include <utilities/ascMalloc.h>
#include <general/tm_time.h>
#include <general/list.h>
#include <solver/slv_types.h>
#include <solver/mtx.h>
#include <solver/var.h>
#include <solver/rel.h>
#include <solver/discrete.h>
#include <solver/conditional.h>
#include <solver/logrel.h>
#include <solver/bnd.h>
#include <solver/slv_common.h>
#include <solver/linsol.h>
#include <solver/linsolqr.h>
#include <solver/linutils.h>
#include <solver/calc.h>
#include <solver/relman.h>
#include <solver/slv_client.h>
#include <solver/system.h>
#include "HelpProc.h"
#include "Sensitivity.h"
#include "HelpProc.h"
#include "SolverGlobals.h"
#include <general/mathmacros.h>

#ifndef lint
static CONST char SensitivityID[] = "$Id: Sensitivity.c,v 1.31 2003/08/23 18:43:08 ballan Exp $";
#endif


#define SAFE_FIX_ME 0

#define DOTIME FALSE
#define DEBUG 1
/*
 * This file attempts to implement the extraction of dy_dx from
 * a system of equations. If one considers a black-box where x are
 * the input variables, u are inuut parameters, y are the output
 * variables, f(x,y,u) is the system of equations that will be solved
 * for given x and u, then one can extract the sensitivity information
 * of y wrt x.
 * One crude but simple way of doing this is to finite-difference the
 * given black box, i.e, perturb x, n\subx times by delta x, resolve
 * the system of equations at the new value of x, and compute
 * dy/dx = (f(x\sub1) - f(x\sub2))/(x\sub1 - x\sub2).
 * This is known to be very expensive.
 *
 * The solution that will be adopted here is to formulate the Jacobian J of
 * the system, (or the system is already solved, to grab the jacobian at
 * the solution. Develop the sensitivity matrix df/dx by exact numnerical
 * differentiation; this will be a n x n\subx matrix. Compute the LU factors
 * for J, and to solve column by column to : LU dz/dx = df/dx. Here
 * z, represents all the internal variables to the system and the strictly
 * output variables y. In other words J = df/dz.
 *
 * Given the solution of the above problem we then simply extract the rows
 * of dz/dx, that pertain to the y variables that we are interested in;
 * this will give us dy/dx.
 */

#ifdef THIS_MAY_BE_UNUSED_CODE
static real64 *zero_vector(real64 *vec, int size)
{
  int c;
  for (c=0;c<size;c++) {
    vec[c] = 0.0;
  }
  return vec;
}
#endif

static int NumberIncludedRels(slv_system_t sys)
{
  static int nrels = -1;
  rel_filter_t rfilter;
  int tmp;

  if (!sys) { /* a NULL system may be used to reinitialise the count */
    nrels = -1;
    return -1;
  }
  if (nrels < 0) {
    /*get used (included) relation count -- equalities only !*/
    rfilter.matchbits = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
    rfilter.matchvalue = (REL_INCLUDED | REL_EQUALITY | REL_ACTIVE);
    tmp = slv_count_solvers_rels(sys,&rfilter);
    nrels = tmp;
    return tmp;
  } else {
    return nrels;
  }
}

static int NumberFreeVars(slv_system_t sys)
{
  static int nvars = -1;
  var_filter_t vfilter;
  int tmp;

  if (!sys) {
    nvars = -1;
    return -1;
  }
  if (nvars < 0) {
    /*get used (free and incident) variable count */

    vfilter.matchbits = (VAR_FIXED |VAR_INCIDENT| VAR_SVAR | VAR_ACTIVE);
    vfilter.matchvalue = (VAR_INCIDENT | VAR_SVAR | VAR_ACTIVE);
    tmp = slv_count_solvers_vars(sys,&vfilter);
    nvars = tmp;
    return tmp;
  } else {
    return nvars;
  }
}


/*
 *********************************************************************
 * The following bit of code does the computation of the matrix
 * dz/dx. It accepts a slv_system_t and a list of 'input' variables.
 * The matrix df/dx now exists and sits to the right of the main
 * square region of the Jacobian. We will do the following in turn
 * for each variable in this list:
 *
 * 1)   Get the variable index, and use it to extract the required column
 *      from the main gradient matrix, this will be stored in a temporary
 *      vector.
 * 2)   We will then clear the column in the original matrix, as we want to
 *  store the caluclated results back in place.
 * 3)   Add the vector extracted in 1) as rhs to the main matrix.
 * 4)   Call linsol solve on this rhs to yield a solution which represents
 *  a column of dx/dx.
 * 6)   Store this vector back in the location cleared out in 2).
 * 7)   Goto 1.
 *
 * At the end of this procedure we would have calculated all the columns of
 * dz/dx, and left them back in the main matrix.
 *********************************************************************
 */


/*
 * At this point we should have an empty jacobian. We now
 * need to call relman_diff over the *entire* matrix.
 * Calculates the entire jacobian. It is initially unscaled.
 *
 * Note: this assumes the sys given is using one of the ascend solvers
 * with either linsol or linsolqr. Both are now allowed. baa 7/95
 */
static int Compute_J(slv_system_t sys)
{
  int32 row;
  var_filter_t vfilter;
  linsol_system_t lin_sys = NULL;
  linsolqr_system_t lqr_sys = NULL;
  mtx_matrix_t mtx;
  struct rel_relation **rlist;
  int nrows;
  int qr=0;
#if DOTIME
  double time1;
#endif

#if DOTIME
  time1 = tm_cpu_time();
#endif
  /*
   * Get the linear system from the solve system.
   * Get the matrix from the linear system.
   * Get the relations list from the solve system.
   */
  lin_sys = slv_get_linsol_sys(sys);
  if (lin_sys==NULL) {
    qr=1;
    lqr_sys=slv_get_linsolqr_sys(sys);
  }
  mtx = slv_get_sys_mtx(sys);
  if (mtx==NULL || (lin_sys==NULL && lqr_sys==NULL)) {
    FPRINTF(stderr,"Compute_dy_dx: error, found NULL.\n");
    return 1;
  }
  rlist = slv_get_solvers_rel_list(sys);
  nrows = NumberIncludedRels(sys);

  calc_ok = TRUE;
  vfilter.matchbits = VAR_SVAR;
  vfilter.matchvalue = VAR_SVAR;
  /*
   * Clear the entire matrix and then compute
   * the gradients at the current point.
   */
  mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
  for(row=0; row<nrows; row++) {
    struct rel_relation *rel;
    /* added  */
    double resid;

    rel = rlist[mtx_row_to_org(mtx,row)];
    (void)relman_diffs(rel,&vfilter,mtx,&resid,SAFE_FIX_ME);

    /* added  */
    rel_set_residual(rel,resid);

  }
  if (qr) {
    linsolqr_matrix_was_changed(lqr_sys);
  } else {
    linsol_matrix_was_changed(lin_sys);
  }
#if DOTIME
  time1 = tm_cpu_time() - time1;
  FPRINTF(stderr,"Time taken for Compute_J = %g\n",time1);
#endif
  return(!calc_ok);
}


/*
 * Note a rhs would have been previously added
 * to keep the system happy.
 * Now handles both linsol and linsolqr as needed.
 * Assumes region to be factored is in upper left corner.
 * Region is reordered using the last reorder method used in
 * the case of linsolqr, so all linsolqr methods are available
 * to this routine.
 */
static int LUFactorJacobian(slv_system_t sys)
{
  linsol_system_t lin_sys=NULL;
  linsolqr_system_t lqr_sys=NULL;
  mtx_region_t region;
  int size,nvars,nrels;
#if DOTIME
  double time1;
#endif

#if DOTIME
  time1 = tm_cpu_time();
#endif

  nvars = NumberFreeVars(sys);
  nrels = NumberIncludedRels(sys);
  size = MAX(nvars,nrels);      /* get size of matrix */
  mtx_region(&region,0,size-1,0,size-1);    /* set the region */
  lin_sys = slv_get_linsol_sys(sys);    /* get the linear system */
  if (lin_sys!=NULL) {
    linsol_matrix_was_changed(lin_sys);
    linsol_reorder(lin_sys,&region);    /* This was entire_MATRIX */
    linsol_invert(lin_sys,&region);     /* invert the given matrix over
                                         * the given region */
  } else {
    enum factor_method fm;
    int oldtiming;
    lqr_sys = slv_get_linsolqr_sys(sys);
    /* WE are ASSUMING that the system has been qr_preped before now! */
    linsolqr_matrix_was_changed(lqr_sys);
    linsolqr_reorder(lqr_sys,&region,natural); /* should retain original */
    fm = linsolqr_fmethod(lqr_sys);
    if (fm == unknown_f) {
      fm = ranki_kw2; /* make sure somebody set it */
    }
    oldtiming = g_linsolqr_timing;
    g_linsolqr_timing = 0;
    linsolqr_factor(lqr_sys,fm);
    g_linsolqr_timing = oldtiming;
  }

#if DOTIME
  time1 = tm_cpu_time() - time1;
  FPRINTF(stderr,"Time taken for LUFactorJacobian = %g\n",time1);
#endif
  return 0;
}


/*
 * At this point the rhs would have already been added.
 *
 * Extended to handle either factorization code:
 * linsol or linsolqr.
 */
static int Compute_dy_dx_smart(slv_system_t sys,
                               real64 *rhs,
                               real64 **dy_dx,
                               int *inputs, int ninputs,
                               int *outputs, int noutputs)
{
  linsol_system_t lin_sys=NULL;
  linsolqr_system_t lqr_sys=NULL;
  mtx_matrix_t mtx;
  int col,current_col;
  int row;
  int capacity;
  real64 *solution = NULL;
  int i,j;
#if DOTIME
  double time1;
#endif

#if DOTIME
  time1 = tm_cpu_time();
#endif

  lin_sys = slv_get_linsol_sys(sys);    /* get the linear system */
  lqr_sys = slv_get_linsolqr_sys(sys);  /* get the linear system */
  mtx = slv_get_sys_mtx(sys);       /* get the matrix */

  capacity = mtx_capacity(mtx);
  solution = (real64 *)asccalloc(capacity,sizeof(real64));

  /*
   * The array inputs is a list of original indexes, of the variables
   * that we are trying to obtain the sensitivity to. We have to
   * get the *current* column from the matrix based on those indices.
   * Hence we use mtx_org_to_col. This little manipulation is not
   * necessary for the computed solution as the solve routine returns
   * the results in the *original* order rather than the *current* order.
   */
  if (lin_sys!=NULL) {
    for (j=0;j<ninputs;j++) {
      col = inputs[j];
      current_col = mtx_org_to_col(mtx,col);
      mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS); /* get the column */

      linsol_rhs_was_changed(lin_sys,rhs);
      linsol_solve(lin_sys,rhs);            /* solve */
      linsol_copy_solution(lin_sys,rhs,solution);   /* get the solution */

      for (i=0;i<noutputs;i++) {    /* copy the solution to dy_dx */
        row = outputs[i];
        dy_dx[i][j] = -1.0*solution[row]; /* the -1 comes from the lin alg */
      }
      /*
       * zero the vector using the incidence pattern of our column.
       * This is faster than the naiive approach of zeroing
       * everything especially if the vector is large.
       */
      mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
    }
  } else {
    for (j=0;j<ninputs;j++) {
      col = inputs[j];
      current_col = mtx_org_to_col(mtx,col);
      mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS);

      linsolqr_rhs_was_changed(lqr_sys,rhs);
      linsolqr_solve(lqr_sys,rhs);
      linsolqr_copy_solution(lqr_sys,rhs,solution);

      for (i=0;i<noutputs;i++) {
        row = outputs[i];
        dy_dx[i][j] = -1.0*solution[row];
      }
      mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
    }
  }
  if (solution) {
    ascfree((char *)solution);
  }

#if DOTIME
  time1 = tm_cpu_time() - time1;
  FPRINTF(stderr,"Time for Compute_dy_dx_smart = %g\n",time1);
#endif
  return 0;
}


/*
 *********************************************************************
 * This code is provided for the benefit of a temporary
 * fix for the derivative problem in Lsode.
 * The proper permanent fix for lsode is to dump it in favor of
 * cvode or dassl.
 * Extended 7/95 baa to deal with linsolqr and blsode.
 * It is assumed the system has been solved at the current point.
 *********************************************************************
 */
int Asc_BLsodeDerivatives(slv_system_t sys, double **dy_dx,
                       int *inputs_ndx_list, int ninputs,
                       int *outputs_ndx_list, int noutputs)
{
  static int n_calls = 0;
  linsolqr_system_t lqr_sys;    /* stuff for the linear system & matrix */
  mtx_matrix_t mtx;
  int32 capacity;
  real64 *scratch_vector = NULL;
  int result=0;

  (void)NumberFreeVars(NULL);       /* used to re-init the system */
  (void)NumberIncludedRels(NULL);   /* used to re-init the system */
  if (!sys) {
    FPRINTF(stderr,"The solve system does not exist !\n");
    return 1;
  }

  result = Compute_J(sys);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
    return 1;
  }

  lqr_sys = slv_get_linsolqr_sys(sys);  /* get the linear system */
  if (lqr_sys==NULL) {
    FPRINTF(stderr,"Early termination due to missing linsolqr system.\n");
    return 1;
  }
  mtx = slv_get_sys_mtx(sys);   /* get the matrix */
  if (mtx==NULL) {
    FPRINTF(stderr,"Early termination due to missing mtx in linsolqr.\n");
    return 1;
  }
  capacity = mtx_capacity(mtx);
  scratch_vector = (real64 *)asccalloc(capacity,sizeof(real64));
  linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);

  result = LUFactorJacobian(sys);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
    goto error;
  }
  result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
                               inputs_ndx_list, ninputs,
                               outputs_ndx_list, noutputs);

  linsolqr_remove_rhs(lqr_sys,scratch_vector);
  if (result) {
    FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
    goto error;
  }

 error:
  n_calls++;
  if (scratch_vector) {
    ascfree((char *)scratch_vector);
  }
  return result;
}


int Asc_MtxNormsCmd(ClientData cdata, Tcl_Interp *interp,
                int argc, CONST84 char *argv[])
{
  slv_system_t sys;
  linsol_system_t lin_sys = NULL;
  linsolqr_system_t linqr_sys = NULL;
  mtx_matrix_t mtx;
  mtx_region_t reg;
  double norm;
  int solver;

  (void)cdata;    /* stop gcc whine about unused parameter */
  (void)argv;     /* stop gcc whine about unused parameter */

  if (argc!=1) {
    Tcl_SetResult(interp, "wrong # args: Usage __mtx_norms", TCL_STATIC);
    return TCL_ERROR;
  }
  sys = g_solvsys_cur;
  if (sys==NULL) {
    Tcl_SetResult(interp, "__mtx_norms called with slv_system", TCL_STATIC);
    return TCL_ERROR;
  }
  solver = slv_get_selected_solver(sys);
  switch(solver) {
  default:
  case 0:
    lin_sys = slv_get_linsol_sys(sys);
    mtx = linsol_get_matrix(lin_sys);
    reg.col.low = reg.row.low = 0;
    reg.col.high = reg.row.high = mtx_symbolic_rank(mtx);
    norm = linutils_A_1_norm(mtx,&reg);
    FPRINTF(stderr,"A_1_norm = %g\n",norm);
    norm = linutils_A_infinity_norm(mtx,&reg);
    FPRINTF(stderr,"A_infinity_norm = %g\n",norm);
    norm = linutils_A_Frobenius_norm(mtx,&reg);
    FPRINTF(stderr,"A_Frobenius_norm = %g\n",norm);
    norm = linutils_A_cond_kaa(lin_sys,mtx,NULL);
    FPRINTF(stderr,"A_condition # = %g\n",norm);
    break;
  case 3:
  case 5:
    linqr_sys = slv_get_linsolqr_sys(sys);
    mtx = linsolqr_get_matrix(linqr_sys);
    reg.col.low = reg.row.low = 0;
    reg.col.high = reg.row.high = mtx_symbolic_rank(mtx);
    norm = linutils_A_1_norm(mtx,&reg);
    FPRINTF(stderr,"A_1_norm = %g\n",norm);
    norm = linutils_A_infinity_norm(mtx,&reg);
    FPRINTF(stderr,"A_infinity_norm = %g\n",norm);
    norm = linutils_A_Frobenius_norm(mtx,&reg);
    FPRINTF(stderr,"A_Frobenius_norm = %g\n",norm);
    norm = linutils_A_condqr_kaa(linqr_sys,mtx,NULL);
    FPRINTF(stderr,"A_condition # = %g\n",norm);
    break;
  }
  return TCL_OK;
}
1	/*
2	* Sensitivity.c
3	* by Kirk Abbott and Ben Allan
4	* Created: 1/94
5	* Version: $Revision: 1.31 $
6	* Version control file: $RCSfile: Sensitivity.c,v $
7	* Date last modified: $Date: 2003/08/23 18:43:08 $
8	* Last modified by: $Author: ballan $
9	*
10	* This file is part of the ASCEND Tcl/Tk interface
11	*
12	* Copyright 1997, Carnegie Mellon University
13	*
14	* The ASCEND Tcl/Tk interface is free software; you can redistribute
15	* it and/or modify it under the terms of the GNU General Public License as
16	* published by the Free Software Foundation; either version 2 of the
17	* License, or (at your option) any later version.
18	*
19	* The ASCEND Tcl/Tk interface is distributed in hope that it will be
20	* useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
21	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22	* General Public License for more details.
23	*
24	* You should have received a copy of the GNU General Public License
25	* along with the program; if not, write to the Free Software Foundation,
26	* Inc., 675 Mass Ave, Cambridge, MA 02139 USA. Check the file named
27	* COPYING. COPYING is found in ../compiler.
28	*/
29
30
31	/*
32	* Sensititvity analysis code.
33	*/
34
35	#include <math.h>
36	#include <tcl.h>
37	#include <utilities/ascConfig.h>
38	#include <utilities/ascMalloc.h>
39	#include <general/tm_time.h>
40	#include <general/list.h>
41	#include <solver/slv_types.h>
42	#include <solver/mtx.h>
43	#include <solver/var.h>
44	#include <solver/rel.h>
45	#include <solver/discrete.h>
46	#include <solver/conditional.h>
47	#include <solver/logrel.h>
48	#include <solver/bnd.h>
49	#include <solver/slv_common.h>
50	#include <solver/linsol.h>
51	#include <solver/linsolqr.h>
52	#include <solver/linutils.h>
53	#include <solver/calc.h>
54	#include <solver/relman.h>
55	#include <solver/slv_client.h>
56	#include <solver/system.h>
57	#include "HelpProc.h"
58	#include "Sensitivity.h"
59	#include "HelpProc.h"
60	#include "SolverGlobals.h"
61	#include <general/mathmacros.h>
62
63	#ifndef lint
64	static CONST char SensitivityID[] = "$Id: Sensitivity.c,v 1.31 2003/08/23 18:43:08 ballan Exp $";
65	#endif
66
67
68	#define SAFE_FIX_ME 0
69
70	#define DOTIME FALSE
71	#define DEBUG 1
72	/*
73	* This file attempts to implement the extraction of dy_dx from
74	* a system of equations. If one considers a black-box where x are
75	* the input variables, u are inuut parameters, y are the output
76	* variables, f(x,y,u) is the system of equations that will be solved
77	* for given x and u, then one can extract the sensitivity information
78	* of y wrt x.
79	* One crude but simple way of doing this is to finite-difference the
80	* given black box, i.e, perturb x, n\subx times by delta x, resolve
81	* the system of equations at the new value of x, and compute
82	* dy/dx = (f(x\sub1) - f(x\sub2))/(x\sub1 - x\sub2).
83	* This is known to be very expensive.
84	*
85	* The solution that will be adopted here is to formulate the Jacobian J of
86	* the system, (or the system is already solved, to grab the jacobian at
87	* the solution. Develop the sensitivity matrix df/dx by exact numnerical
88	* differentiation; this will be a n x n\subx matrix. Compute the LU factors
89	* for J, and to solve column by column to : LU dz/dx = df/dx. Here
90	* z, represents all the internal variables to the system and the strictly
91	* output variables y. In other words J = df/dz.
92	*
93	* Given the solution of the above problem we then simply extract the rows
94	* of dz/dx, that pertain to the y variables that we are interested in;
95	* this will give us dy/dx.
96	*/
97
98	#ifdef THIS_MAY_BE_UNUSED_CODE
99	static real64 zero_vector(real64 vec, int size)
100	{
101	int c;
102	for (c=0;c<size;c++) {
103	vec[c] = 0.0;
104	}
105	return vec;
106	}
107	#endif
108
109	static int NumberIncludedRels(slv_system_t sys)
110	{
111	static int nrels = -1;
112	rel_filter_t rfilter;
113	int tmp;
114
115	if (!sys) { /* a NULL system may be used to reinitialise the count */
116	nrels = -1;
117	return -1;
118	}
119	if (nrels < 0) {
120	/get used (included) relation count -- equalities only !/
121	rfilter.matchbits = (REL_INCLUDED \| REL_EQUALITY \| REL_ACTIVE);
122	rfilter.matchvalue = (REL_INCLUDED \| REL_EQUALITY \| REL_ACTIVE);
123	tmp = slv_count_solvers_rels(sys,&rfilter);
124	nrels = tmp;
125	return tmp;
126	} else {
127	return nrels;
128	}
129	}
130
131	static int NumberFreeVars(slv_system_t sys)
132	{
133	static int nvars = -1;
134	var_filter_t vfilter;
135	int tmp;
136
137	if (!sys) {
138	nvars = -1;
139	return -1;
140	}
141	if (nvars < 0) {
142	/get used (free and incident) variable count /
143
144	vfilter.matchbits = (VAR_FIXED \|VAR_INCIDENT\| VAR_SVAR \| VAR_ACTIVE);
145	vfilter.matchvalue = (VAR_INCIDENT \| VAR_SVAR \| VAR_ACTIVE);
146	tmp = slv_count_solvers_vars(sys,&vfilter);
147	nvars = tmp;
148	return tmp;
149	} else {
150	return nvars;
151	}
152	}
153
154
155	/*
156	*********************************************************************
157	* The following bit of code does the computation of the matrix
158	* dz/dx. It accepts a slv_system_t and a list of 'input' variables.
159	* The matrix df/dx now exists and sits to the right of the main
160	* square region of the Jacobian. We will do the following in turn
161	* for each variable in this list:
162	*
163	* 1) Get the variable index, and use it to extract the required column
164	* from the main gradient matrix, this will be stored in a temporary
165	* vector.
166	* 2) We will then clear the column in the original matrix, as we want to
167	* store the caluclated results back in place.
168	* 3) Add the vector extracted in 1) as rhs to the main matrix.
169	* 4) Call linsol solve on this rhs to yield a solution which represents
170	* a column of dx/dx.
171	* 6) Store this vector back in the location cleared out in 2).
172	* 7) Goto 1.
173	*
174	* At the end of this procedure we would have calculated all the columns of
175	* dz/dx, and left them back in the main matrix.
176	*********************************************************************
177	*/
178
179
180	/*
181	* At this point we should have an empty jacobian. We now
182	* need to call relman_diff over the entire matrix.
183	* Calculates the entire jacobian. It is initially unscaled.
184	*
185	* Note: this assumes the sys given is using one of the ascend solvers
186	* with either linsol or linsolqr. Both are now allowed. baa 7/95
187	*/
188	static int Compute_J(slv_system_t sys)
189	{
190	int32 row;
191	var_filter_t vfilter;
192	linsol_system_t lin_sys = NULL;
193	linsolqr_system_t lqr_sys = NULL;
194	mtx_matrix_t mtx;
195	struct rel_relation **rlist;
196	int nrows;
197	int qr=0;
198	#if DOTIME
199	double time1;
200	#endif
201
202	#if DOTIME
203	time1 = tm_cpu_time();
204	#endif
205	/*
206	* Get the linear system from the solve system.
207	* Get the matrix from the linear system.
208	* Get the relations list from the solve system.
209	*/
210	lin_sys = slv_get_linsol_sys(sys);
211	if (lin_sys==NULL) {
212	qr=1;
213	lqr_sys=slv_get_linsolqr_sys(sys);
214	}
215	mtx = slv_get_sys_mtx(sys);
216	if (mtx==NULL \|\| (lin_sys==NULL && lqr_sys==NULL)) {
217	FPRINTF(stderr,"Compute_dy_dx: error, found NULL.\n");
218	return 1;
219	}
220	rlist = slv_get_solvers_rel_list(sys);
221	nrows = NumberIncludedRels(sys);
222
223	calc_ok = TRUE;
224	vfilter.matchbits = VAR_SVAR;
225	vfilter.matchvalue = VAR_SVAR;
226	/*
227	* Clear the entire matrix and then compute
228	* the gradients at the current point.
229	*/
230	mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
231	for(row=0; row<nrows; row++) {
232	struct rel_relation *rel;
233	/* added */
234	double resid;
235
236	rel = rlist[mtx_row_to_org(mtx,row)];
237	(void)relman_diffs(rel,&vfilter,mtx,&resid,SAFE_FIX_ME);
238
239	/* added */
240	rel_set_residual(rel,resid);
241
242	}
243	if (qr) {
244	linsolqr_matrix_was_changed(lqr_sys);
245	} else {
246	linsol_matrix_was_changed(lin_sys);
247	}
248	#if DOTIME
249	time1 = tm_cpu_time() - time1;
250	FPRINTF(stderr,"Time taken for Compute_J = %g\n",time1);
251	#endif
252	return(!calc_ok);
253	}
254
255
256	/*
257	* Note a rhs would have been previously added
258	* to keep the system happy.
259	* Now handles both linsol and linsolqr as needed.
260	* Assumes region to be factored is in upper left corner.
261	* Region is reordered using the last reorder method used in
262	* the case of linsolqr, so all linsolqr methods are available
263	* to this routine.
264	*/
265	static int LUFactorJacobian(slv_system_t sys)
266	{
267	linsol_system_t lin_sys=NULL;
268	linsolqr_system_t lqr_sys=NULL;
269	mtx_region_t region;
270	int size,nvars,nrels;
271	#if DOTIME
272	double time1;
273	#endif
274
275	#if DOTIME
276	time1 = tm_cpu_time();
277	#endif
278
279	nvars = NumberFreeVars(sys);
280	nrels = NumberIncludedRels(sys);
281	size = MAX(nvars,nrels); /* get size of matrix */
282	mtx_region(&region,0,size-1,0,size-1); /* set the region */
283	lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
284	if (lin_sys!=NULL) {
285	linsol_matrix_was_changed(lin_sys);
286	linsol_reorder(lin_sys,&region); /* This was entire_MATRIX */
287	linsol_invert(lin_sys,&region); /* invert the given matrix over
288	* the given region */
289	} else {
290	enum factor_method fm;
291	int oldtiming;
292	lqr_sys = slv_get_linsolqr_sys(sys);
293	/* WE are ASSUMING that the system has been qr_preped before now! */
294	linsolqr_matrix_was_changed(lqr_sys);
295	linsolqr_reorder(lqr_sys,&region,natural); /* should retain original */
296	fm = linsolqr_fmethod(lqr_sys);
297	if (fm == unknown_f) {
298	fm = ranki_kw2; /* make sure somebody set it */
299	}
300	oldtiming = g_linsolqr_timing;
301	g_linsolqr_timing = 0;
302	linsolqr_factor(lqr_sys,fm);
303	g_linsolqr_timing = oldtiming;
304	}
305
306	#if DOTIME
307	time1 = tm_cpu_time() - time1;
308	FPRINTF(stderr,"Time taken for LUFactorJacobian = %g\n",time1);
309	#endif
310	return 0;
311	}
312
313
314	/*
315	* At this point the rhs would have already been added.
316	*
317	* Extended to handle either factorization code:
318	* linsol or linsolqr.
319	*/
320	static int Compute_dy_dx_smart(slv_system_t sys,
321	real64 *rhs,
322	real64 **dy_dx,
323	int *inputs, int ninputs,
324	int *outputs, int noutputs)
325	{
326	linsol_system_t lin_sys=NULL;
327	linsolqr_system_t lqr_sys=NULL;
328	mtx_matrix_t mtx;
329	int col,current_col;
330	int row;
331	int capacity;
332	real64 *solution = NULL;
333	int i,j;
334	#if DOTIME
335	double time1;
336	#endif
337
338	#if DOTIME
339	time1 = tm_cpu_time();
340	#endif
341
342	lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
343	lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
344	mtx = slv_get_sys_mtx(sys); /* get the matrix */
345
346	capacity = mtx_capacity(mtx);
347	solution = (real64 *)asccalloc(capacity,sizeof(real64));
348
349	/*
350	* The array inputs is a list of original indexes, of the variables
351	* that we are trying to obtain the sensitivity to. We have to
352	* get the current column from the matrix based on those indices.
353	* Hence we use mtx_org_to_col. This little manipulation is not
354	* necessary for the computed solution as the solve routine returns
355	* the results in the original order rather than the current order.
356	*/
357	if (lin_sys!=NULL) {
358	for (j=0;j<ninputs;j++) {
359	col = inputs[j];
360	current_col = mtx_org_to_col(mtx,col);
361	mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS); /* get the column */
362
363	linsol_rhs_was_changed(lin_sys,rhs);
364	linsol_solve(lin_sys,rhs); /* solve */
365	linsol_copy_solution(lin_sys,rhs,solution); /* get the solution */
366
367	for (i=0;i<noutputs;i++) { /* copy the solution to dy_dx */
368	row = outputs[i];
369	dy_dx[i][j] = -1.0solution[row]; / the -1 comes from the lin alg */
370	}
371	/*
372	* zero the vector using the incidence pattern of our column.
373	* This is faster than the naiive approach of zeroing
374	* everything especially if the vector is large.
375	*/
376	mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
377	}
378	} else {
379	for (j=0;j<ninputs;j++) {
380	col = inputs[j];
381	current_col = mtx_org_to_col(mtx,col);
382	mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS);
383
384	linsolqr_rhs_was_changed(lqr_sys,rhs);
385	linsolqr_solve(lqr_sys,rhs);
386	linsolqr_copy_solution(lqr_sys,rhs,solution);
387
388	for (i=0;i<noutputs;i++) {
389	row = outputs[i];
390	dy_dx[i][j] = -1.0*solution[row];
391	}
392	mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
393	}
394	}
395	if (solution) {
396	ascfree((char *)solution);
397	}
398
399	#if DOTIME
400	time1 = tm_cpu_time() - time1;
401	FPRINTF(stderr,"Time for Compute_dy_dx_smart = %g\n",time1);
402	#endif
403	return 0;
404	}
405
406
407	/*
408	*********************************************************************
409	* This code is provided for the benefit of a temporary
410	* fix for the derivative problem in Lsode.
411	* The proper permanent fix for lsode is to dump it in favor of
412	* cvode or dassl.
413	* Extended 7/95 baa to deal with linsolqr and blsode.
414	* It is assumed the system has been solved at the current point.
415	*********************************************************************
416	*/
417	int Asc_BLsodeDerivatives(slv_system_t sys, double **dy_dx,
418	int *inputs_ndx_list, int ninputs,
419	int *outputs_ndx_list, int noutputs)
420	{
421	static int n_calls = 0;
422	linsolqr_system_t lqr_sys; /* stuff for the linear system & matrix */
423	mtx_matrix_t mtx;
424	int32 capacity;
425	real64 *scratch_vector = NULL;
426	int result=0;
427
428	(void)NumberFreeVars(NULL); /* used to re-init the system */
429	(void)NumberIncludedRels(NULL); /* used to re-init the system */
430	if (!sys) {
431	FPRINTF(stderr,"The solve system does not exist !\n");
432	return 1;
433	}
434
435	result = Compute_J(sys);
436	if (result) {
437	FPRINTF(stderr,"Early termination due to failure in calc Jacobian\n");
438	return 1;
439	}
440
441	lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
442	if (lqr_sys==NULL) {
443	FPRINTF(stderr,"Early termination due to missing linsolqr system.\n");
444	return 1;
445	}
446	mtx = slv_get_sys_mtx(sys); /* get the matrix */
447	if (mtx==NULL) {
448	FPRINTF(stderr,"Early termination due to missing mtx in linsolqr.\n");
449	return 1;
450	}
451	capacity = mtx_capacity(mtx);
452	scratch_vector = (real64 *)asccalloc(capacity,sizeof(real64));
453	linsolqr_add_rhs(lqr_sys,scratch_vector,FALSE);
454
455	result = LUFactorJacobian(sys);
456	if (result) {
457	FPRINTF(stderr,"Early termination due to failure in LUFactorJacobian\n");
458	goto error;
459	}
460	result = Compute_dy_dx_smart(sys, scratch_vector, dy_dx,
461	inputs_ndx_list, ninputs,
462	outputs_ndx_list, noutputs);
463
464	linsolqr_remove_rhs(lqr_sys,scratch_vector);
465	if (result) {
466	FPRINTF(stderr,"Early termination due to failure in Compute_dy_dx\n");
467	goto error;
468	}
469
470	error:
471	n_calls++;
472	if (scratch_vector) {
473	ascfree((char *)scratch_vector);
474	}
475	return result;
476	}
477
478
479	int Asc_MtxNormsCmd(ClientData cdata, Tcl_Interp *interp,
480	int argc, CONST84 char *argv[])
481	{
482	slv_system_t sys;
483	linsol_system_t lin_sys = NULL;
484	linsolqr_system_t linqr_sys = NULL;
485	mtx_matrix_t mtx;
486	mtx_region_t reg;
487	double norm;
488	int solver;
489
490	(void)cdata; /* stop gcc whine about unused parameter */
491	(void)argv; /* stop gcc whine about unused parameter */
492
493	if (argc!=1) {
494	Tcl_SetResult(interp, "wrong # args: Usage __mtx_norms", TCL_STATIC);
495	return TCL_ERROR;
496	}
497	sys = g_solvsys_cur;
498	if (sys==NULL) {
499	Tcl_SetResult(interp, "__mtx_norms called with slv_system", TCL_STATIC);
500	return TCL_ERROR;
501	}
502	solver = slv_get_selected_solver(sys);
503	switch(solver) {
504	default:
505	case 0:
506	lin_sys = slv_get_linsol_sys(sys);
507	mtx = linsol_get_matrix(lin_sys);
508	reg.col.low = reg.row.low = 0;
509	reg.col.high = reg.row.high = mtx_symbolic_rank(mtx);
510	norm = linutils_A_1_norm(mtx,&reg);
511	FPRINTF(stderr,"A_1_norm = %g\n",norm);
512	norm = linutils_A_infinity_norm(mtx,&reg);
513	FPRINTF(stderr,"A_infinity_norm = %g\n",norm);
514	norm = linutils_A_Frobenius_norm(mtx,&reg);
515	FPRINTF(stderr,"A_Frobenius_norm = %g\n",norm);
516	norm = linutils_A_cond_kaa(lin_sys,mtx,NULL);
517	FPRINTF(stderr,"A_condition # = %g\n",norm);
518	break;
519	case 3:
520	case 5:
521	linqr_sys = slv_get_linsolqr_sys(sys);
522	mtx = linsolqr_get_matrix(linqr_sys);
523	reg.col.low = reg.row.low = 0;
524	reg.col.high = reg.row.high = mtx_symbolic_rank(mtx);
525	norm = linutils_A_1_norm(mtx,&reg);
526	FPRINTF(stderr,"A_1_norm = %g\n",norm);
527	norm = linutils_A_infinity_norm(mtx,&reg);
528	FPRINTF(stderr,"A_infinity_norm = %g\n",norm);
529	norm = linutils_A_Frobenius_norm(mtx,&reg);
530	FPRINTF(stderr,"A_Frobenius_norm = %g\n",norm);
531	norm = linutils_A_condqr_kaa(linqr_sys,mtx,NULL);
532	FPRINTF(stderr,"A_condition # = %g\n",norm);
533	break;
534	}
535	return TCL_OK;
536	}