From maintainers-request at octave dot org Wed Nov 17 15:42:02 2004
Subject: Re: Octave 2.1.63 available for ftp
From: David Bateman <David dot Bateman at motorola dot com>
To: maintainers at octave dot org
Date: Wed, 17 Nov 2004 22:41:44 +0100

> In addition to these changes, I think we still have the following
> items for Octave 3.0:
>
>  * allow [x{:}] = f (...) to work
>
>  * private functions
>
>  * some support for MPI
>
>  * bring the manual up to date

If we are discussing the wish-list for 3.0 again, I'd suggest that if
private functions are in the classes should also be added, since the
two are closely related. Classes are just private directories for a
certain data-type. Once such a private direcory for a data-type
existing adding the overload functions won't be too difficult.

Another thing you might need in that case is the dispatch code from
octave-forge, rewritten to fit into the class structure. In any case,
unless you intend to add the sparse code to all of the functions like
lu, inv, chol, qr, det, ranks, sum, reshape, atan2, diag, etc, etc,
etc, then a proper dispatch/class structure will be needed.

Another thing that I'd wish for, for version 3, would be that the
FIXES and patches directories of octave-forge are gone through and the
stuff there either updated and ported to the version 3.0 or definitely
dropped. These directories shouldn't exist in octave-forge.

It would also be nice if the "\" and "/" operators could recognize a
triangular matrix and call the appropriate lapack/blas code. The stuff in
octave-forge for this is a bit of a kludge as it only works for chol and
doesn't use the lapack code for "\" and "/". Maybe something like

   enum {diagonal = 0, upper, lower, full} matrix_type;
   matrix_type typ = diagonal;

   int i, r, c;
   for (i = 0; i < numel (); i++)
     if (elem (i) != 0.)
       {
	 r = i % nr;
   	 c = i / nc;
	 if (r != c)
	   break;
       }

   if (i != numel ())
     if (r < c)
       {
         type = lower;
         for (int j = i + 1; j < numel (); j++)
           {
	     int r2 = i2 % nr;
   	     int c2 = i2 / nc;
	     if (r2 > c2)
	       {
	         type = full;
	         break;
	       }
	   }
       }
     else
       {
         type = upper;
         for (int j = i + 1; j < numel (); j++)
           {
	     int r2 = i2 % nr;
   	     int c2 = i2 / nc;
	     if (r2 < c2)
	      {
	        type = full;
	        break;
	      }
	  }
       }
   }

   if (type == lower || type == lower)
     {
       char typ = 'L';
       if (type = upper)
         type = 'U';

       retval = Matrix (b);
       double *result = retval.fortran_vec ();
       const double *a_data = a.fortran_vec ();
       double one = 1.0;
       char job = '1';
       F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
                                  F77_CONST_CHAR_ARG2 (&typ, 1),
                                  F77_CONST_CHAR_ARG2 ("N", 1),
				  nc, tmp_data, nr, 
				  rcond, pz, piz, info
				  F77_CHAR_ARG_LEN (1)));
	      
       if (f77_exception_encountered)
	(*current_liboctave_error_handler) 
	  ("unrecoverable error in dtrcon");
	      
       if (info != 0) 
	 info = -2;

       volatile double rcond_plus_one = rcond + 1.0;

       if (rcond_plus_one == 1.0 || xisnan (rcond))
	 {
	   info = -2;

	  if (sing_handler)
	    sing_handler (rcond);
	  else
	    (*current_liboctave_error_handler)
	      ("matrix singular to machine precision, rcond = %g",
	       rcond);
	 }
       else
	 {
           F77_XFCN (dtrsm, DTRSM, (F77_CONST_CHAR_ARG2 ("L", 1),
                     F77_CONST_CHAR_ARG2 (&typ, 1),
                     F77_CONST_CHAR_ARG2 ("N", 1),
                     F77_CONST_CHAR_ARG2 ("N", 1), nr, b_nc, 
                     one, a_data, nr, result, b.rows ()));

           if (f77_exception_encountered)
	     (*current_liboctave_error_handler) ("unrecoverable error in dgetrf")
         }
     }
   else
     {
       // existing solve code

       ...
     }

this might be added to the Matrix::solve and ComplexMatrix::solve
functions. I don't think this would be very costly, as in most cases
you'd quickly identify full matrices. It would be interest to see how
this performs ralative to the chol function in octave-forge.

I also think version 3.0 should be the release where octave-forge is 
decruftified; there is lots of stuff to support older releases, that
should be dropped at that point, though that has nothing to do with
octave per-se...

Cheers
David

-- 
David Bateman                                David dot Bateman at motorola dot com
Motorola CRM                                 +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary