where is ipvt (from base-lu) defined?

[John W. Eaton]

On 29-Mar-2005, Clinton Chee <address@hidden> wrote:

| Hi Octave Maintainers,
| 
| I'm trying to track down a problem with initialization of ipvt in base-lu.cc
| 
| Version: Octave 2.1.64
| 
| Platform: Itanium2
| Linux 2.4.21-sgi303r2 #1 SMP Mon Nov 29 15:29:38 PST 2004 ia64 ia64 ia64
| GNU/Linux
| 
| Compiler:
| Reading specs from /usr/lib/gcc-lib/ia64-redhat-linux/3.2.3/specs
| Configured with: ../configure --prefix=/usr --mandir=/usr/share/man
| --infodir=/usr/share/info --enable-shared --enable-threads=posix
| --disable-checking --with-system-zlib --enable-__cxa_atexit
| --host=ia64-redhat-linux
| Thread model: posix
| gcc version 3.2.3 20030502 (Red Hat Linux 3.2.3-34)
| 
| NOTE: The problem occurs on the above platform but DOES NOT occur on a
| 32bit Linux Cluster.
| 
| Description:
| 1. Start octave
| 2. type: a=([1,2;3,4])
| 3. type: lu(a)
| 
| The program segfaults, the  gdb backtrace is given below. I narrowed the
| problem down to base-lu.cc in the "P" function. In the "P" function,
| there is:
| 
| int k = ipvt.xelem(i)
| 
| .... within a for loop where i=0,1.
| 
| On the Itanium2 platform,
| ipvt.xelem(0)=1
| ipvt.xelem(1)=-1
| 
| On the 32 bit linux cluster
| ipvt.xelem(0)=1
| ipvt.xelem(1)=1
| 
| I understand why it Seg Faults (k given the value of -1, and then
| pvt.xelem(k) fails)
| 
| 
| Question: Where does ipvt become defined (in specific, when you do a
| "lu(a)" operation)? I know it is declared in base-lu.h as an MArray<T>
| type. But why on one platforms it has values of {1,1} and on another {1,-1}

It's defined in the LU constructor.  For example, in dbleLU.cc:

LU::LU (const Matrix& a)
{
  int a_nr = a.rows ();
  int a_nc = a.cols ();
  int mn = (a_nr < a_nc ? a_nr : a_nc);

  ipvt.resize (mn);
  int *pipvt = ipvt.fortran_vec ();

  a_fact = a;
  double *tmp_data = a_fact.fortran_vec ();

  int info = 0;

  F77_XFCN (dgetrf, DGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));

  if (f77_exception_encountered)
    (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
  else
    ipvt -= 1;
}

First it's resized, then values are set for it in the call to DGETRF.
If the calculation succeeds, the values in ipvt are decremented so
that they match the zero-based indexing used for arrays in Octave.

jwe