use LAPACK routine for triangular systems?

classic Classic list List threaded Threaded
13 messages Options
Reply | Threaded
Open this post in threaded view
|

use LAPACK routine for triangular systems?

Evan Monroig
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
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

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;
}
Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Paul Kienzle
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.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Bill Denney
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.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

William Poetra Yoga H-2
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.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Evan Monroig
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.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Evan Monroig
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.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

John W. Eaton-6
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.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Evan Monroig-2
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.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Stéfan van der Walt
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?DaCodaAlFine

The wiki also has information on wrapping Fortran code at

http://wiki.octave.org/wiki.pl?OctaveFortran

Regards
Stéfan



-------------------------------------------------------------
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
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Fredrik Lingvall
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
Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Paul Kienzle
In reply to this post by Evan Monroig

On Nov 27, 2005, at 8:22 PM, Evan Monroig wrote:

> 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 ;)

In octave-forge/extra/linear-algebra there is a triangular
matrix type and implements a solver.

You will have to check out the 2.9.x sparse support yourself,
but I believe David implemented optimizations for banded
and perhaps triangular matrices.

- Paul



-------------------------------------------------------------
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
-------------------------------------------------------------

Reply | Threaded
Open this post in threaded view
|

Re: use LAPACK routine for triangular systems?

Evan Monroig
On 11/28/05, Paul Kienzle <[hidden email]> wrote:

>
> On Nov 27, 2005, at 8:22 PM, Evan Monroig wrote:
>
> > 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 ;)
>
> In octave-forge/extra/linear-algebra there is a triangular
> matrix type and implements a solver.

Yes, I just looked at it and saw that
octave-forge/extra/linear-algebra/ov-re-tri.cc implements a triangular
matrix type, and also a method for solving a linear system with such a
matrix by recursive back-substitution (I'm not sure if this is the
right name).

After reading the corresponding functions in LAPACK and BLAS
(lapack3-3.0.20000531a/blas/src/dtrsm.f), the same method seems to be
used. I don't quite fluently read fortran but it looked like that.

I thought that there would be special methods to avoid troubles with
badly-conditionned matrices when dividing by the diagonal numbers of
the matrix ^^.

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
-------------------------------------------------------------