help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Paul Kienzle
Subject: Re: A(i,:) in .oct file?
Date: Fri, 9 Feb 2001 07:55:42 +0000
User-agent: Mutt/1.2.5i

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, and the code
runs 2-3 times faster if I implement the BLAS directly in C++.  I've
included daxpy and dscal below.  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?

Paul Kienzle
address@hidden

On Thu, Feb 08, 2001 at 10:34:31PM -0600, A S Hodel wrote:
> You could possibly use a direct call the the BLAS daxpy (or zaxpy) routine.
> I think it would go something like:
> 
>   daxpy( vector_len, -l(i), g(i,0), g.rows(), g(i+1,0), g.rows() );
> 
> That depends on Octave storing its matrices in column-major order, but it's
> the best I can come up with.
> 
> > on 2/8/01 9:37 PM, John W. Eaton at address@hidden wrote:
> 
> > On  7-Feb-2001, Paul Kienzle <address@hidden> wrote:
> > 
> > | Hi!  How do I translate something like:
> > | 
> > |  g(i+1,:) -= l(i)*g(i,:);
> > | 
> > | into a .oct file?
> > | 
> > |  for (j=0; j<g.columns(); j++) g(i+1,j) -= l(i)*g(i,j);
> > | 
> > | works, but it isn't very elegant,
> > 
> > I can't think of a better way.  Assuming g is a Matrix object, I think
> > you could use
> > 
> > g.insert (g.row(i+1) - l(i)*g.row(i), i+1, 0);
> > 
> > but that is going to generate some copies.  Probably the loop is about
> > as good as you can do.
> > 
> > | and I do have a lot of them to do.
> > 
> > I suppose I would define a macro or inline function to hide the
> > details.
> > 
> > jwe
> > 
> 
> -- 
> A S Hodel Assoc. Prof. Dept Elec. Comp. Eng.  Auburn Univ, AL 36849-5201
> (334) 844-1854 Fax -1809 http://www.eng.auburn.edu/~scotte
> 
> 
------
Some C++blas.  I grant these to the public domain.

// daxpy: y <- ax + y
// daxpy(n,a,&x,dx,&y,dy)
//    n>0 is the length of the vector
//    a!=0 is the constant multiplier
//    x is the source vector
//    dx is the source increment
//    y is the update vector
//    dy is the update increment
void 
daxpy (const int n, const double a, const double x[], 
       const int dx, double y[], const int dy)
{
  if (n > 0 && a != 0)
    if (dx == 1 && dy == 1)
      {
        const int m = n%4;
        for (int i=0; i < m; i++) y[i] += a*x[i];
        for (int i=m; i < n; i+=4) 
          {
            y[i] += a*x[i];
            y[i+1] += a*x[i+1];
            y[i+2] += a*x[i+2];
            y[i+3] += a*x[i+3];
          }
      }
    else
      {
        int ix = (dx < 0 ? (-n+1)*dx : 0);
        int iy = (dy < 0 ? (-n+1)*dy : 0);
        for (int i=0; i < n; i++, ix+=dx, iy+=dy) y[iy] += a*x[ix];
      }
}


// dscal: y <- ay
// dscal(n,a,&y,dy)
//    n>0 is the length of the vector
//    a is the constant multiplier
//    y is the update vector
//    dy>0 is the update increment
void 
dscal (const int n, const double a, double y[], const int dy)
{
  if (n > 0 && dy > 0)
    if (dy == 1)
      {
        const int m = n%5;
        for (int i=0; i < m; i++) y[i] *= a;
        for (int i=m; i < n; i+=5) 
          {
            y[i] *= a;
            y[i+1] *= a;
            y[i+2] *= a;
            y[i+3] *= a;
            y[i+4] *= a;
          }
      }
    else
      {
        const int last=n*dy;
        for (int i=0; i < last; i+=dy) y[i] *= a;
      }
}



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



reply via email to

[Prev in Thread] Current Thread [Next in Thread]