Re: use LAPACK routine for triangular systems?

[Fredrik Lingvall]

Evan Monroig wrote:

> Hi,
>
> I know that octave uses LAPACK routines for some linear functions (QR
> factorization, etc).
>
> In my algorithm I end with a triangular system to solve, and I know
> that LAPACK has some specific routines for triangular systems (like
> this one http://www.netlib.org/lapack/double/dtrtrs.f). Is there a way
> to take advantage of these with octave?
>
> Evan
>
>
>  
I have made I few mex-files for some of the BLAS and LAPACK routines
for that purpose. I have attached  an example  using the LAPACK POTRI
routine. You can use the mex-tools in octave-forge to build the 
corresponding
oct-files for octave (you need to link to your BLAS and LAPACK libs).

/Fredrik

#include "mex.h"
#include <math.h>
#include <string.h>

/***
 *
 *  LAPACK Positive definite invert using factorization (POTRI).
 *
 *
 * Fredrik Lingvall 2003-06-12  
 *
 ***/

//
// Function prototypes.
//



// DPOTRF computes the Cholesky factorization of a real symmetric positive 
definite matrix A.
extern void dpotrf_(const char *UPLO, const int *N, 
                    const double *A, const int *lda,
                    const int *info);

extern void dpotri_(const char *UPLO, const int *N, 
                    const double *A, const int *lda,
                    const int *info);


/*** From dpotri.f ******
     

SUBROUTINE DPOTRI( UPLO, N, A, LDA, INFO )
*
*  -- LAPACK routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     March 31, 1993
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INFO, LDA, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * )
*     ..
*
*  Purpose
*  =======
*
*  DPOTRI computes the inverse of a real symmetric positive definite
*  matrix A using the Cholesky factorization A = U**T*U or A = L*L**T
*  computed by DPOTRF.
*
*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          = 'U':  Upper triangle of A is stored;


*          = 'U':  Upper triangle of A is stored;
*          = 'L':  Lower triangle of A is stored.
*
*  N       (input) INTEGER
*          The order of the matrix A.  N >= 0.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
*          On entry, the triangular factor U or L from the Cholesky
*          factorization A = U**T*U or A = L*L**T, as computed by
*          DPOTRF.
*          On exit, the upper or lower triangle of the (symmetric)
*          inverse of A, overwriting the input factor U or L.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, the (i,i) element of the factor U or L is
*                zero, and the inverse could not be computed.
*
*  =====================================================================
*/


/***
 * 
 * Matlab (MEX) gateway function
 *
 ***/


void  mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  // Input arguments.
  int          M, N, lda, info;
  register int nn, kk, K;
  double       *A, *Y;
  char         *upper_or_lower = "L";
  
  //  
  // Check for proper number of arguments
  //

  if (nrhs != 1) {
    mexErrMsgTxt("One input argument required!");
  }

  if (nlhs > 1) {
    mexErrMsgTxt("Too many output arguments, one expected!");
  }
  
  //
  // Check array geometry args.
  //

  M = mxGetM(prhs[0]);
  N = mxGetN(prhs[0]);
  A = mxGetPr(prhs[0]);

  if (mxIsSparse(prhs[0]))
    mexErrMsgTxt("Input marix cannot be sparse!");

  if (M != N)
    mexErrMsgTxt("Marix must be square!");
  
  // Allocate mamory for output matrix.
  plhs[0] = mxCreateDoubleMatrix(M,N, mxREAL);
  Y = mxGetPr(plhs[0]);
  if (!Y)
    mexErrMsgTxt("Memory allocation failed!");
  
  // Copy A matrix (since Y is overwritten by DPOTRI).
  memcpy(Y, A, M*N*sizeof(double));

  lda = N;                      // Leading dim in A.
  dpotrf_(upper_or_lower, &N, Y, &lda, &info);
  dpotri_(upper_or_lower, &N, Y, &lda, &info);

  //printf("info = %d\n",info);

  // Copy lower triangle to upper triangle (keep iteration variables in CPU 
regs).
  K = M;
  for (kk=1; kk<K; kk++)
    for (nn=0; nn<kk; nn++)
      Y[nn+kk*K] = Y[kk+nn*K];

  return;
}