[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: use LAPACK routine for triangular systems?
From: |
Fredrik Lingvall |
Subject: |
Re: use LAPACK routine for triangular systems? |
Date: |
Fri, 25 Nov 2005 10:08:44 +0100 |
User-agent: |
Mozilla Thunderbird 1.0.7 (X11/20051019) |
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;
}
Re: use LAPACK routine for triangular systems?, Paul Kienzle, 2005/11/25