# use LAPACK routine for triangular systems?

13 messages
Open this post in threaded view
|

## use LAPACK routine for triangular systems?

 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 ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 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 #include /***  *  *  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
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 In reply to this post by Evan Monroig Octave-forge has trisolve.m in main/splines which uses the LAPACK routines.  Octave 2.9.x has this built into the sparse matrix solver routines. - Paul On Nov 24, 2005, at 2:33 AM, 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 > > > > ------------------------------------------------------------- > Octave is freely available under the terms of the GNU GPL. > > Octave's home on the web:  http://www.octave.org> How to fund new projects:  http://www.octave.org/funding.html> Subscription information:  http://www.octave.org/archive.html> ------------------------------------------------------------- > ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 I don't know how it is currently done, but may we be able to generate code from netlib code automatically into octave code?  I think that I recall (though feel free to correct me if I'm wrong) that most netlib code is generated automatically from templates?  If it is, how difficult would it be for someone in the octave community to add automatic octave code generation to those templates? Bill On Fri, 25 Nov 2005, Paul Kienzle wrote: > Octave-forge has trisolve.m in main/splines which uses the LAPACK > routines.  Octave 2.9.x has this built into the sparse matrix > solver routines. > > - Paul > > On Nov 24, 2005, at 2:33 AM, 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 >> >> >> >> ------------------------------------------------------------- >> Octave is freely available under the terms of the GNU GPL. >> >> Octave's home on the web:  http://www.octave.org>> How to fund new projects:  http://www.octave.org/funding.html>> Subscription information:  http://www.octave.org/archive.html>> ------------------------------------------------------------- >> > > > > ------------------------------------------------------------- > Octave is freely available under the terms of the GNU GPL. > > Octave's home on the web:  http://www.octave.org> How to fund new projects:  http://www.octave.org/funding.html> Subscription information:  http://www.octave.org/archive.html> ------------------------------------------------------------- > -- "I am not one of those who believe that a great army is the means of maintaining peace, because if you build up a great profession those who form parts of it want to exercise their profession."    -- Woodrow Wilson ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 On 11/26/05, Bill Denney <[hidden email]> wrote: > I don't know how it is currently done, but may we be able to generate code > from netlib code automatically into octave code?  I think that I recall > (though feel free to correct me if I'm wrong) that most netlib code is > generated automatically from templates?  If it is, how difficult would it > be for someone in the octave community to add automatic octave code > generation to those templates? > > Bill > If I recall correctly, Octave already includes some netlib code (that are public domain). Look in the subdirectory libcruft in the Octave sources. -- William Poetra Yoga Hadisoeseno ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 In reply to this post by Fredrik Lingvall On 11/25/05, Fredrik Lingvall <[hidden email]> wrote: > >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? > > > 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). Thanks, I will try your method. If I understand well, the process to do that is: 1) get the LAPACK f-file 2) wrap up the file in the way you did to have a mex-file 3) use mex-tools to convert that to an oct-file Am I right? Evan ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 In reply to this post by Paul Kienzle On 11/26/05, Paul Kienzle <[hidden email]> wrote: > Octave-forge has trisolve.m in main/splines which uses the LAPACK > routines.  Octave 2.9.x has this built into the sparse matrix > solver routines. > > - Paul This is for tridiagonal systems, and I am looking into using a LAPACK routine for upper-triangular systems ;) Thanks anyway ! ^^ Evan ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 In reply to this post by Evan Monroig On 28-Nov-2005, Evan Monroig wrote: | 1) get the LAPACK f-file | 2) wrap up the file in the way you did to have a mex-file | 3) use mex-tools to convert that to an oct-file If you are starting from scratch, write the code required to make a .oct file directly instead of using Matlab's mex interface. jwe ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 On 11/28/05, John W. Eaton <[hidden email]> wrote: > On 28-Nov-2005, Evan Monroig wrote: > > | 1) get the LAPACK f-file > | 2) wrap up the file in the way you did to have a mex-file > | 3) use mex-tools to convert that to an oct-file > > If you are starting from scratch, write the code required to make a > .oct file directly instead of using Matlab's mex interface. Ok, thanks ! After a little googling I found this link: http://octave.sourceforge.net/coda/Is this the "best-practice" to do the conversion? Evan ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 On Mon, Nov 28, 2005 at 02:29:14PM +0900, Evan Monroig wrote: > Ok, thanks ! > > After a little googling I found this link: > http://octave.sourceforge.net/coda/The preferred version of Da Coda Al Fine is on the Octave Wiki at http://wiki.octave.org/wiki.pl?DaCodaAlFineThe wiki also has information on wrapping Fortran code at http://wiki.octave.org/wiki.pl?OctaveFortranRegards Stéfan ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web:  http://www.octave.orgHow to fund new projects:  http://www.octave.org/funding.htmlSubscription information:  http://www.octave.org/archive.html-------------------------------------------------------------
Open this post in threaded view
|

## Re: use LAPACK routine for triangular systems?

 In reply to this post by Evan Monroig Evan Monroig wrote: ```On 11/25/05, Fredrik Lingvall [hidden email] wrote: ``` ```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? ``` ```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). ``` ``` Thanks, I will try your method. If I understand well, the process to do that is: 1) get the LAPACK f-file ``` You really don't need the LAPACK f-file. You link the LAPACK that probably is on your system already. /F