Re: 2.9.12 (2.9.13) chol broken with --enable-64?

[Fredrik Lingvall]

John W. Eaton wrote:
> On 11-Aug-2007, Fredrik Lingvall wrote:
>
> | I have compiled LAPACK (3.0.1) with -m64 and -march=nocona and I use 
>
> What compiler?  
That was given in the original post, but again:

fllap5 src # gcc -v
Using built-in specs.
Target: x86_64-pc-linux-gnu
Configured with: 
/var/tmp/portage/sys-devel/gcc-4.1.2/work/gcc-4.1.2/configure 
--prefix=/usr --bindir=/usr/x86_64-pc-linux-gnu/gcc-bin/4.1.2 
--includedir=/usr/lib/gcc/x86_64-pc-linux-gnu/4.1.2/include 
--datadir=/usr/share/gcc-data/x86_64-pc-linux-gnu/4.1.2 
--mandir=/usr/share/gcc-data/x86_64-pc-linux-gnu/4.1.2/man 
--infodir=/usr/share/gcc-data/x86_64-pc-linux-gnu/4.1.2/info 
--with-gxx-include-dir=/usr/lib/gcc/x86_64-pc-linux-gnu/4.1.2/include/g++-v4 
--host=x86_64-pc-linux-gnu --build=x86_64-pc-linux-gnu --disable-altivec 
--enable-nls --without-included-gettext --with-system-zlib 
--disable-checking --disable-werror --enable-secureplt 
--disable-libunwind-exceptions --enable-multilib --enable-libmudflap 
--disable-libssp --disable-libgcj --enable-languages=c,c++,fortran 
--enable-shared --enable-threads=posix --enable-__cxa_atexit 
--enable-clocale=gnu
Thread model: posix
gcc version 4.1.2 (Gentoo 4.1.2 p1.0.1)

> What do those options mean?
>
> | GOTO BLAS with
> | BINARY64  = 1 in the Makefile.rule file.
>
> I don't know what that does.
>   
I'm not 100% sure either, but GOTO's lib is rather common so maybe 
someone else on this list nows?
> Octave expects to use the Fortran interface to lapack and blas
> functions, so do these options make the INTEGER arguments to those
> functions signed 8 byte values?
>   
From gcc's man page:

`-m32'
`-m64'
    Generate code for a 32-bit or 64-bit environment.  The 32-bit
    environment sets int, long and pointer to 32 bits and generates
    code that runs on any i386 system.  The 64-bit environment sets
    int to 32 bits and long and pointer to 64 bits and generates code
    for AMD's x86-64 architecture.

> | To test I created an oct file that calls the LAPACK routines,
> | 
> |   dpotrf_(upper_or_lower, &N, Y, &lda, &info);  (cholesky factorization)
> | 
> | and  
> | 
> |   dpotri_(upper_or_lower, &N, Y, &lda, &info); (positive triangualar 
> | inverse)
> | 
> | which should do the same thing as cholinv (I guess). The result is
> | 
> | octave:1> A=rand(100,100);
> | octave:2> A=A'*A;
> | octave:3> C = potri(A);
> | octave:4> C = potri(A);
> | octave:5> C = potri(A);
>
> What is potri? 
>   
That is the oct-file that contains the calls to dpotrf_(...) and  
dpotri_(...), respectively. Its a simple wrapper
to the LAPACK Fortran rountines  above (I have attached the file now).
> | octave:6> cholinv(A);
> | octave:7> cholinv(A);
> | error: cholinv: matrix not positive definite
> | 
> | I test the info parameter after both calls above and the info parameter 
> | is always = 0 (no error).
>
> I don't think anyone can help you with this problem if you don't show
> the code you are using.
>   
From what I can see (in the liboctave/dbleCHOL.cc), octave's chol and 
cholinv also uses dpotrf_(....) and 
dpotri_(...).  I have assumed,

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

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

that is, N and lda are int's. Perhaps, they should be long's instead? 
Or, maybe the gfortan
option -fdefault-integer-8 should be used when compiling LAPACK?

/Fredrik


#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>

//
// Octave headers.
//

#include <octave/oct.h>

#include <octave/config.h>

#include <iostream>
using namespace std;

#include <octave/defun-dld.h>
#include <octave/error.h>
#include <octave/oct-obj.h>
#include <octave/pager.h>
#include <octave/symtab.h>
#include <octave/variables.h>

/***
 *
 *  LAPACK Positive definite invert using factorization (POTRI).
 *
 *
 * Fredrik Lingvall 2007-08-11 : File created.  
 *
 ***/

//
// Function prototypes.
//



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

extern "C" 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.
*
*  =====================================================================
*/

/***
 * 
 * Octave (Oct) gateway function for dgemm.
 *
 ***/

DEFUN_DLD (potri, args, nlhs,
           "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {}  @var{y} = potri (@var{a},@var{s}).\n\
\n\
potri : Positive definite invert using factorization.\n\
\n\
@var{a}     - An NxN positive definite symmetric matrix.\n\
@var{s}     - 'L' (lower triangular) or 'U' (upper triangular). Optional 
parameter.\n\
\n\
Output parameter:\n\
\n\
@var{y} - Output matrix (optional).\n\
\n\
potri is an oct-function.\n\
\n\
@copyright{ 2006-08-25 Fredrik Lingvall}.\n\
@end deftypefn\n\
@seealso {BLAS and LAPACK documentation}")
{
  int    M, N, lda, info, VERBOSE=0, len, in_place=0;
  int    n, k, K, buflen;
  double *A, *Y;
  char   upper_or_lower[50];

  octave_value_list oct_retval; 
  
  int nrhs = args.length ();
  
  //  
  // Check for proper number of arguments
  //
  
  if (nrhs > 2) {
    error("Too many input arguments!");
    return oct_retval;
  }

  if (nlhs > 1) {
    error("Too many output arguments, one expected!");
    return oct_retval;
  }

  if (nrhs == 2) {

    //if (!mxIsChar(prhs[1])) {
    if  (!args(1).is_string()) {
      error("2nd argument must be a string (L or U) ");
      return oct_retval;
    }
    
    //len = mxGetN(prhs[1]) + mxGetM(prhs[1])+1;

    // Copy the 2-byte string.
    //mxGetString(prhs[1], upper_or_lower, len);
    std::string strin = args(1).string_value();     
    buflen = strin.length();
    for ( n=0; n<buflen; n++ ) {
      upper_or_lower[n] = strin[n];
    }
    

    if ( strcmp(upper_or_lower,"L") == 0 ) {
      in_place = 0;
      if (VERBOSE)
        printf("Using lower triangular computation.\n");

    } else if ( strcmp(upper_or_lower,"Li") == 0 ) {
      in_place = 1;
      if (nlhs == 1) {
        error("No output argument expected for in-place operation!");
        return oct_retval;
      }
      if (VERBOSE)
        printf("Using lower triangular in-place computation.\n");

    } else if ( strcmp(upper_or_lower,"U") == 0) {
      in_place = 0;
      if (VERBOSE)
        printf("Using upper triangular computation.\n");

    } else if ( strcmp(upper_or_lower,"Ui") == 0) {
      in_place = 1;
      if (nlhs == 1) {
        error("No output argument expected for in-place operation!");
        return oct_retval;
      }
      if (VERBOSE)
        printf("Using upper triangular in-place computation.\n");

    } else {
      error("Unknown string in 2nd  argument! Must be L (lower) or U (upper).");
      return oct_retval;
    }
    
  } else { // Use lower triangular by default.
    strcpy(upper_or_lower,"L");
    
    if (VERBOSE)
      printf("Using lower triangular computation.\n");
  }


  //
  // Check array geometry args.
  //

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

  const Matrix tmp = args(0).matrix_value();
  M = tmp.rows();
  N = tmp.cols();
  A = (double*) tmp.fortran_vec();

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

  if (M != N) {
    error("Marix must be square!");
    return oct_retval;
  }

  if (in_place) { // In-place operation.
    Y = A;
  } else {
    // Allocate mamory for output matrix.
    //plhs[0] = mxCreateDoubleMatrix(M,N, mxREAL);
    //Y = mxGetPr(plhs[0]);
    Matrix Ymat(M, N);
    Y = Ymat.fortran_vec();
    oct_retval.append(Ymat);

    if (!Y) {
      error("Memory allocation failed!");
      return oct_retval;
    }
    // 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);

  if (info !=0)
    printf("info = %d\n",info);
  
  dpotri_(upper_or_lower, &N, Y, &lda, &info);

  if (info !=0)
    printf("info = %d\n",info);

  if (nrhs == 1) { // Use the full matrix.
    // Copy lower triangle to upper triangle.
    K = M;
    for (k=1; k<K; k++)
      for (n=0; n<k; n++)
        Y[n+k*K] = Y[k+n*K];
  }
  return oct_retval;
}