/************************  Using C routines in Gauss under OS/2

Author: Richard Stanton
Date: 3 Sep 1996
Provided without guarantees for public non-commercial use.


I recently spent a while working out how to compile C code into a
DLL file using the emx/gcc compiler under OS/2, and how to use the
result with GAUSS. Doing this for just one small routine resulted in a
7-fold speedup in my program.

1) Here is the C code itself, rhs.c.  It contains just one routine,
tribksub1, which I use for repeatedly solving the equation Ax = d, for
different values of the vector d, and where A is tridiagonal. For
those who are interested, I have previously calculated the LU
decomposition of A, and this routine now performs forward/backward
substitution to calculate x. Since each value affects the next, it has
to be done as a loop, which is why it's so slow in GAUSS.
***************************************************************/


/* rhs.c (emx+gcc) DLL for GAUSS */

int tribksub1(double n[], double a[], double beta[], double gamma[], double d[],
                double x[])
/*
 * Solve equation A x = d using output of routine TRILU (beta, gamma).
 * Suited for multiple RHS, since no divisions are performed in this routine.
 */
{
        int i;
/*
 * Forward substitution
 */
        x[0] = d[0] * beta[0];
        for (i=1; i<n[0]; i++) {
                x[i] = (d[i] - a[i] * x[i-1]) * beta[i];
        }
/*
 * Backward substitution
 */
        for (i=n[0]-2; i>=0; i--) {
                x[i] -= gamma[i+1] * x[i+1];
        }
        return 0;
}



/*****************************************************************
2) Here's the file rhs.def, which tells the compiler what function
   names to export, and some other stuff which I don't know or care
   about. It works, so here it is:
******************************************************************/


;
; rhs.def
;
LIBRARY RHS INITINSTANCE TERMINSTANCE
DESCRIPTION 'rhs'
DATA
   MULTIPLE NONSHARED
EXPORTS
   tribksub1    @1



/*************************************************************
3) To create the file rhs.dll, use the command

   gcc -Zdll -Zomf -Zso -Zsys rhs.c rhs.def

4) Here's the section of GAUSS code that calls the function:

   dlibrary e:\gauss\dlib\rhs.dll;
   dllcall tribksub1(ny1, a, lu1, lu2, dL, prices);

   The output is stored in the vector prices.


There are a couple of things to be aware of:

1) As Ron S pointed out, Gauss now has built-in functions to solve
   equations of the type I was solving, which would probably be even
   faster.

2) My function doesn't take the actual LU
   decomposition, but at least one of the vectors has all its elements
   inverted, or something similar, so that when I do the substitutions
   to solve the equation, I only use multiplication, and no division
   (which is typically a slower operation).

The above provides an example of how to use DLLs with GAUSS under
OS/2. If you actually want it to be usable, here is the
function trilu, which calculates the (modified) LU decomposition used
as input to tribksub1.

Richard Stanton
*************************************************************************/



/*
 * Calculate modified LU decomposition of tridiagonal matrix
 * as input for TRIBKSUB, solution of tridiagonal simultaneous equ.
 * Key to parameters passed:
 * a: array of left of diagonal entries. (Must include a[1] = .0)
 * b: array of diagonal entries.
 * c: array of right of diagonal entries. (Must include c[n] = .0)
 *
 * Call as {lu1, lu2} = trilu(a,b,c)
 *
 * lu1, lu2: outputs of routine, used with a as subsequent input for
 * TRIBKSUB.
 */
proc (2) = trilu(a, b, c);

  local i, n, lu1, lu2, temp;

  if b[1] == .0;
    errorlog "Error 1 in TRILU";
    end;
  endif;

  n = rows(a);

  lu1 = zeros(n,1);
  lu2 = zeros(n,1);

  lu1[1] = 1.0 / b[1];
  lu2[1] = 0.0;
  i = 2;
  do while i <= n;
    lu2[i] = c[i-1] * lu1[i-1];
    temp = b[i] - a[i] * lu2[i];
    if temp == 0.0;
      errorlog "Error 2 in TRILU";
      end;
    endif;
    lu1[i] = 1.0 / temp;
    i = i+1;
  endo;

  retp(lu1, lu2);
endp;