Re: A(i,:) in .oct file?

[John W. Eaton]

On  9-Feb-2001, Paul Kienzle <address@hidden> wrote:

| Thanks for the suggestion.  I was able to call the fortran blas directly
| using:
| 
|       #include <octave/f77-fcn.h>
|       extern "C" 
|       {
|         void F77_FCN (daxpy, DAXPY) 
|           (const int& n, const double& a, const double x[], 
|            const int& dx, double y[], const int& dy);
|       }
| 
|       inline void 
|       daxpy (const int n, const double a, const double x[], 
|              const int dx, double y[], const int dy)
|       {
|         F77_XFCN (daxpy, DAXPY, (n, a, x, dx, y, dy));
|       }
| 
| Doug, you might want to provide this example in your oct-file FAQ for
| anyone who wants to include Fortran code in their projects.  Note
| particularly that all fortran arguments are passed by reference!

| However, the Fortran calling interface has some overhead,

You can lower that to almost nothing by using F77_FCN instead of
F77_XFCN.  The latter macro is only needed if Fortran code you are
calling can ever call XSTOPX, which in turn calls various things that
eventually result in a longjmp call back to the place where F77_XFCN
was used.  Daxpy doesn't need this level of protection.

| These can be used for the m-file constructs
| 
|       g(i,:) += a*g(j,:); 
|       g(i,:) *= a;
| 
| as follows:
| 
|       const int n = g.rows();
|       const int c = g.columns();
|       daxpy(c, a, &g(j,0), n, &g(i,0), n);
|       dscal(c, a, &g(i,0));
| 
| Something more for the oct-file FAQ?

Please be sure to note that this could modify more than one Matrix
object, and that you are probably better off doing something like

  double *gptr = g.fortran_vec ();

and then doing the offset calculations yourself (the call to
fortran_vec forces the actual data handled by g to not be shared).
Also, note that using &g(j,0) and passing n as the increment and
expecting daxpy to work on a row of g assumes column-major order for
the storage, which may not always be true if the design of the classes
changes in the future.

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