generic/packages/sensitivity.c

/*  ASCEND modelling environment
    Copyright (C) 1990-2006 Carnegie-Mellon University

    This program 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, or (at your option)
    any later version.

    This program is distributed in the 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 this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330,
    Boston, MA 02111-1307, USA.
*//** @file
    Sensititvity analysis code

    @see documentation in sensitivity.h.
*//*
    by Kirk Abbott
*/

#include "sensitivity.h"

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <utilities/error.h>
#include <utilities/ascMalloc.h>
#include <utilities/ascPanic.h>
#include <compiler/packages.h>
#include <compiler/fractions.h>
#include <compiler/dimen.h>
#include <compiler/atomvalue.h>
#include <compiler/instquery.h>
#include <general/mathmacros.h>

#define DEBUG 1

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

/*----------------------------------------------------
  UTILITY ROUTINES
*/

/**
    Allocate memory for a matrix
    @param nrows Number of rows
    @param ncols Number of colums
    @return Pointer to the allocated matrix memory location
*/
real64 **make_matrix(int nrows, int ncols){
  real64 **result;
  int i;
  result = (real64 **)calloc(nrows,sizeof(real64*));
  for (i=0;i<nrows;i++) {
    result[i] = (real64 *)calloc(ncols,sizeof(real64));
  }
  return result;
}

/**
    Free a matrix from memory
    @param matrix Memory location for the matrix
    @param nrows Number of rows in the matrix
*/
void free_matrix(real64 **matrix, int nrows){
  int i;
  if (!matrix)
    return;
  for (i=0;i<nrows;i++) {
    if (matrix[i]) {
      free(matrix[i]);
      matrix[i] = NULL;
    }
  }
  free(matrix);
}

/**
    Fetch an element from a branch of an arglist...?
*/
struct Instance *FetchElement(struct gl_list_t *arglist,
                     unsigned long whichbranch,
                     unsigned long whichelement){
  struct gl_list_t *branch;
  struct Instance *element;

  if (!arglist) return NULL;
  branch = (struct gl_list_t *)gl_fetch(arglist,whichbranch);
  element = (struct Instance *)gl_fetch(branch,whichelement);
  return element;
}

/**
    Looks like it returns the number of relations in a solver's system

    @param sys The system in question
    @return int Number of relations in the system
    @see slv_count_solvers_rels
*/
int NumberRels(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;
}

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;
  }
}

/**
    Looks like it returns the number of free variables in a solver's system

    @param sys The system in question
    @return the number of free variables

    @see slv_count_solvers_vars
*/
int NumberFreeVars(slv_system_t sys){
  static int nvars = -1;
  var_filter_t vfilter;
  int tmp;

  if (!sys) {
    /* no system: return error */
    nvars = -1;
    return -1;
  }

  if (nvars >=0){
    /* return stored value */
    return nvars;
  }

  /* 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;
}

/*
 * 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.
 * fixed and free.
 *
 *  Calculates the entire jacobian.
 *  It is initially unscaled.
 */
int Compute_J(slv_system_t sys)
{
  int32 row;
  var_filter_t vfilter;
  linsolqr_system_t lqr_sys;
  mtx_matrix_t mtx;
  struct rel_relation **rlist;
  int nrows;
  real64 resid;

  /*
   * Get the linear system from the solve system.
   * Get the matrix from the linear system.
   * Get the relations list from the solve system.
   */
  lqr_sys = slv_get_linsolqr_sys(sys);
  mtx = linsolqr_get_matrix(lqr_sys);
  rlist = slv_get_solvers_rel_list(sys);
  nrows = slv_get_num_solvers_rels(sys);

  calc_ok = TRUE;
  vfilter.matchbits =  (VAR_SVAR | VAR_ACTIVE) ;
  vfilter.matchvalue =  vfilter.matchbits ;

  /*
   * 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;
    rel = rlist[mtx_row_to_org(mtx,row)];
    asc_assert(rel!=NULL);
    (void)relman_diffs(rel,&vfilter,mtx,&resid,1);
  }
  linsolqr_matrix_was_changed(lqr_sys);

  return(!calc_ok);
}

/**
    @NOTE there is another (probably better?) version of this in models/sensitivity
    that uses sparse matrices. This one got pulled out of some files in tcltk dir
    and is used by the LSODE integrator.

 * 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.
*/
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;
  CONSOLE_DEBUG("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.

    This routine is part of the 'temporary solution' for derivatives in lsode.c
*/
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 = ASC_NEW_ARRAY_CLEAR(real64,capacity);

  /*
   * 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;
  CONSOLE_DEBUG("Time for Compute_dy_dx_smart = %g\n",time1);
#endif
  return 0;
}

#undef DEBUG
1	/* ASCEND modelling environment
2	Copyright (C) 1990-2006 Carnegie-Mellon University
3
4	This program is free software; you can redistribute it and/or modify
5	it under the terms of the GNU General Public License as published by
6	the Free Software Foundation; either version 2, or (at your option)
7	any later version.
8
9	This program is distributed in the hope that it will be useful,
10	but WITHOUT ANY WARRANTY; without even the implied warranty of
11	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12	GNU General Public License for more details.
13
14	You should have received a copy of the GNU General Public License
15	along with this program; if not, write to the Free Software
16	Foundation, Inc., 59 Temple Place - Suite 330,
17	Boston, MA 02111-1307, USA.
18	//* @file
19	Sensititvity analysis code
20
21	@see documentation in sensitivity.h.
22	//
23	by Kirk Abbott
24	*/
25
26	#include "sensitivity.h"
27
28	#include <stdio.h>
29	#include <string.h>
30	#include <math.h>
31	#include <utilities/error.h>
32	#include <utilities/ascMalloc.h>
33	#include <utilities/ascPanic.h>
34	#include <compiler/packages.h>
35	#include <compiler/fractions.h>
36	#include <compiler/dimen.h>
37	#include <compiler/atomvalue.h>
38	#include <compiler/instquery.h>
39	#include <general/mathmacros.h>
40
41	#define DEBUG 1
42
43	#if 0
44	static real64 zero_vector(real64 vec, int size)
45	{
46	int c;
47	for (c=0;c<size;c++) {
48	vec[c] = 0.0;
49	}
50	return vec;
51	}
52	#endif
53
54	/*----------------------------------------------------
55	UTILITY ROUTINES
56	*/
57
58	/**
59	Allocate memory for a matrix
60	@param nrows Number of rows
61	@param ncols Number of colums
62	@return Pointer to the allocated matrix memory location
63	*/
64	real64 **make_matrix(int nrows, int ncols){
65	real64 **result;
66	int i;
67	result = (real64 *)calloc(nrows,sizeof(real64));
68	for (i=0;i<nrows;i++) {
69	result[i] = (real64 *)calloc(ncols,sizeof(real64));
70	}
71	return result;
72	}
73
74	/**
75	Free a matrix from memory
76	@param matrix Memory location for the matrix
77	@param nrows Number of rows in the matrix
78	*/
79	void free_matrix(real64 **matrix, int nrows){
80	int i;
81	if (!matrix)
82	return;
83	for (i=0;i<nrows;i++) {
84	if (matrix[i]) {
85	free(matrix[i]);
86	matrix[i] = NULL;
87	}
88	}
89	free(matrix);
90	}
91
92	/**
93	Fetch an element from a branch of an arglist...?
94	*/
95	struct Instance FetchElement(struct gl_list_t arglist,
96	unsigned long whichbranch,
97	unsigned long whichelement){
98	struct gl_list_t *branch;
99	struct Instance *element;
100
101	if (!arglist) return NULL;
102	branch = (struct gl_list_t *)gl_fetch(arglist,whichbranch);
103	element = (struct Instance *)gl_fetch(branch,whichelement);
104	return element;
105	}
106
107	/**
108	Looks like it returns the number of relations in a solver's system
109
110	@param sys The system in question
111	@return int Number of relations in the system
112	@see slv_count_solvers_rels
113	*/
114	int NumberRels(slv_system_t sys){
115	static int nrels = -1;
116	rel_filter_t rfilter;
117	int tmp;
118
119	if (!sys) { /* a NULL system may be used to reinitialise the count */
120	nrels = -1;
121	return -1;
122	}
123	if (nrels < 0) {
124	/get used (included) relation count -- equalities only !/
125	rfilter.matchbits = (REL_INCLUDED \| REL_EQUALITY \| REL_ACTIVE);
126	rfilter.matchvalue = (REL_INCLUDED \| REL_EQUALITY \| REL_ACTIVE);
127	tmp = slv_count_solvers_rels(sys,&rfilter);
128	nrels = tmp;
129	return tmp;
130	}
131	else return nrels;
132	}
133
134	int NumberIncludedRels(slv_system_t sys){
135	static int nrels = -1;
136	rel_filter_t rfilter;
137	int tmp;
138
139	if (!sys) { /* a NULL system may be used to reinitialise the count */
140	nrels = -1;
141	return -1;
142	}
143	if (nrels < 0) {
144	/get used (included) relation count -- equalities only !/
145	rfilter.matchbits = (REL_INCLUDED \| REL_EQUALITY \| REL_ACTIVE);
146	rfilter.matchvalue = (REL_INCLUDED \| REL_EQUALITY \| REL_ACTIVE);
147	tmp = slv_count_solvers_rels(sys,&rfilter);
148	nrels = tmp;
149	return tmp;
150	} else {
151	return nrels;
152	}
153	}
154
155	/**
156	Looks like it returns the number of free variables in a solver's system
157
158	@param sys The system in question
159	@return the number of free variables
160
161	@see slv_count_solvers_vars
162	*/
163	int NumberFreeVars(slv_system_t sys){
164	static int nvars = -1;
165	var_filter_t vfilter;
166	int tmp;
167
168	if (!sys) {
169	/* no system: return error */
170	nvars = -1;
171	return -1;
172	}
173
174	if (nvars >=0){
175	/* return stored value */
176	return nvars;
177	}
178
179	/* get used (free and incident) variable count */
180	vfilter.matchbits = (VAR_FIXED \| VAR_INCIDENT \| VAR_SVAR \| VAR_ACTIVE);
181	vfilter.matchvalue = (VAR_INCIDENT \| VAR_SVAR \| VAR_ACTIVE);
182	tmp = slv_count_solvers_vars(sys,&vfilter);
183	nvars = tmp;
184	return tmp;
185	}
186
187	/*
188	* The following bit of code does the computation of the matrix
189	* dz/dx. It accepts a slv_system_t and a list of 'input' variables.
190	* The matrix df/dx now exists and sits to the right of the main
191	* square region of the Jacobian. We will do the following in turn
192	* for each variable in this list:
193	*
194	* 1) Get the variable index, and use it to extract the required column
195	* from the main gradient matrix, this will be stored in a temporary
196	* vector.
197	* 2) We will then clear the column in the original matrix, as we want to
198	* store the caluclated results back in place.
199	* 3) Add the vector extracted in 1) as rhs to the main matrix.
200	* 4) Call linsol solve on this rhs to yield a solution which represents
201	* a column of dx/dx.
202	* 6) Store this vector back in the location cleared out in 2).
203	* 7) Goto 1.
204	*
205	* At the end of this procedure we would have calculated all the columns of
206	* dz/dx, and left them back in the main matrix.
207	*/
208
209
210	/*
211	* At this point we should have an empty jacobian. We now
212	* need to call relman_diff over the entire matrix.
213	* fixed and free.
214	*
215	* Calculates the entire jacobian.
216	* It is initially unscaled.
217	*/
218	int Compute_J(slv_system_t sys)
219	{
220	int32 row;
221	var_filter_t vfilter;
222	linsolqr_system_t lqr_sys;
223	mtx_matrix_t mtx;
224	struct rel_relation **rlist;
225	int nrows;
226	real64 resid;
227
228	/*
229	* Get the linear system from the solve system.
230	* Get the matrix from the linear system.
231	* Get the relations list from the solve system.
232	*/
233	lqr_sys = slv_get_linsolqr_sys(sys);
234	mtx = linsolqr_get_matrix(lqr_sys);
235	rlist = slv_get_solvers_rel_list(sys);
236	nrows = slv_get_num_solvers_rels(sys);
237
238	calc_ok = TRUE;
239	vfilter.matchbits = (VAR_SVAR \| VAR_ACTIVE) ;
240	vfilter.matchvalue = vfilter.matchbits ;
241
242	/*
243	* Clear the entire matrix and then compute
244	* the gradients at the current point.
245	*/
246	mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
247	for(row=0; row<nrows; row++) {
248	struct rel_relation *rel;
249	rel = rlist[mtx_row_to_org(mtx,row)];
250	asc_assert(rel!=NULL);
251	(void)relman_diffs(rel,&vfilter,mtx,&resid,1);
252	}
253	linsolqr_matrix_was_changed(lqr_sys);
254
255	return(!calc_ok);
256	}
257
258	/**
259	@NOTE there is another (probably better?) version of this in models/sensitivity
260	that uses sparse matrices. This one got pulled out of some files in tcltk dir
261	and is used by the LSODE integrator.
262
263	* Note a rhs would have been previously added
264	* to keep the system happy.
265	* Now handles both linsol and linsolqr as needed.
266	* Assumes region to be factored is in upper left corner.
267	* Region is reordered using the last reorder method used in
268	* the case of linsolqr, so all linsolqr methods are available
269	* to this routine.
270	*/
271	int LUFactorJacobian(slv_system_t sys){
272	linsol_system_t lin_sys=NULL;
273	linsolqr_system_t lqr_sys=NULL;
274	mtx_region_t region;
275	int size,nvars,nrels;
276	#if DOTIME
277	double time1;
278	#endif
279
280	#if DOTIME
281	time1 = tm_cpu_time();
282	#endif
283
284	nvars = NumberFreeVars(sys);
285	nrels = NumberIncludedRels(sys);
286	size = MAX(nvars,nrels); /* get size of matrix */
287	mtx_region(&region,0,size-1,0,size-1); /* set the region */
288	lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
289	if (lin_sys!=NULL) {
290	linsol_matrix_was_changed(lin_sys);
291	linsol_reorder(lin_sys,&region); /* This was entire_MATRIX */
292	linsol_invert(lin_sys,&region); /* invert the given matrix over
293	* the given region */
294	} else {
295	enum factor_method fm;
296	int oldtiming;
297	lqr_sys = slv_get_linsolqr_sys(sys);
298	/* WE are ASSUMING that the system has been qr_preped before now! */
299	linsolqr_matrix_was_changed(lqr_sys);
300	linsolqr_reorder(lqr_sys,&region,natural); /* should retain original */
301	fm = linsolqr_fmethod(lqr_sys);
302	if (fm == unknown_f) {
303	fm = ranki_kw2; /* make sure somebody set it */
304	}
305	oldtiming = g_linsolqr_timing;
306	g_linsolqr_timing = 0;
307	linsolqr_factor(lqr_sys,fm);
308	g_linsolqr_timing = oldtiming;
309	}
310
311	#if DOTIME
312	time1 = tm_cpu_time() - time1;
313	CONSOLE_DEBUG("Time taken for LUFactorJacobian = %g\n",time1);
314	#endif
315	return 0;
316	}
317
318
319	/**
320	At this point the rhs would have already been added.
321
322	Extended to handle either factorization code:
323	linsol or linsolqr.
324
325	This routine is part of the 'temporary solution' for derivatives in lsode.c
326	*/
327	int Compute_dy_dx_smart(slv_system_t sys,
328	real64 *rhs,
329	real64 **dy_dx,
330	int *inputs, int ninputs,
331	int *outputs, int noutputs)
332	{
333	linsol_system_t lin_sys=NULL;
334	linsolqr_system_t lqr_sys=NULL;
335	mtx_matrix_t mtx;
336	int col,current_col;
337	int row;
338	int capacity;
339	real64 *solution = NULL;
340	int i,j;
341	#if DOTIME
342	double time1;
343	#endif
344
345	#if DOTIME
346	time1 = tm_cpu_time();
347	#endif
348
349	lin_sys = slv_get_linsol_sys(sys); /* get the linear system */
350	lqr_sys = slv_get_linsolqr_sys(sys); /* get the linear system */
351	mtx = slv_get_sys_mtx(sys); /* get the matrix */
352
353	capacity = mtx_capacity(mtx);
354	solution = ASC_NEW_ARRAY_CLEAR(real64,capacity);
355
356	/*
357	* The array inputs is a list of original indexes, of the variables
358	* that we are trying to obtain the sensitivity to. We have to
359	* get the current column from the matrix based on those indices.
360	* Hence we use mtx_org_to_col. This little manipulation is not
361	* necessary for the computed solution as the solve routine returns
362	* the results in the original order rather than the current order.
363	*/
364	if (lin_sys!=NULL) {
365	for (j=0;j<ninputs;j++) {
366	col = inputs[j];
367	current_col = mtx_org_to_col(mtx,col);
368	mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS); /* get the column */
369
370	linsol_rhs_was_changed(lin_sys,rhs);
371	linsol_solve(lin_sys,rhs); /* solve */
372	linsol_copy_solution(lin_sys,rhs,solution); /* get the solution */
373
374	for (i=0;i<noutputs;i++) { /* copy the solution to dy_dx */
375	row = outputs[i];
376	dy_dx[i][j] = -1.0solution[row]; / the -1 comes from the lin alg */
377	}
378	/*
379	* zero the vector using the incidence pattern of our column.
380	* This is faster than the naiive approach of zeroing
381	* everything especially if the vector is large.
382	*/
383	mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
384	}
385	} else {
386	for (j=0;j<ninputs;j++) {
387	col = inputs[j];
388	current_col = mtx_org_to_col(mtx,col);
389	mtx_org_col_vec(mtx,current_col,rhs,mtx_ALL_ROWS);
390
391	linsolqr_rhs_was_changed(lqr_sys,rhs);
392	linsolqr_solve(lqr_sys,rhs);
393	linsolqr_copy_solution(lqr_sys,rhs,solution);
394
395	for (i=0;i<noutputs;i++) {
396	row = outputs[i];
397	dy_dx[i][j] = -1.0*solution[row];
398	}
399	mtx_zr_org_vec_using_col(mtx,current_col,rhs,mtx_ALL_ROWS);
400	}
401	}
402	if (solution) {
403	ascfree((char *)solution);
404	}
405
406	#if DOTIME
407	time1 = tm_cpu_time() - time1;
408	CONSOLE_DEBUG("Time for Compute_dy_dx_smart = %g\n",time1);
409	#endif
410	return 0;
411	}
412
413	#undef DEBUG