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;

}