[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: A(i,:) in .oct file?
From: |
John W. Eaton |
Subject: |
Re: A(i,:) in .oct file? |
Date: |
Fri, 9 Feb 2001 12:38:31 -0600 |
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
-------------------------------------------------------------