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 "interface/HelpProc.h"
#include "interface/Sensitivity.h"
#include "interface/HelpProc.h"
#include "interface/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	aw0a	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 "interface/HelpProc.h"
58			#include "interface/Sensitivity.h"
59			#include "interface/HelpProc.h"
60			#include "interface/SolverGlobals.h"
61	johnpye	122	#include "general/mathmacros.h"
62	aw0a	1
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			}