Re: Fortran

[Stef Pillaert]

Op dinsdag 15 mei 2001 16:50, schreef Przemek Klosowski:
>    How can I dinamically link my FORTRAN functions with Octave? In the
> manual only a C++ function is described and I don't know how to write a C++
> wrapper funtion (because I don't know C++....).
Here 's asomewhat complex example: I call the dsysvx - function from the 
fortran lapack-library (to solve symmetric systems) My system is 
Redhat-linux. I'm not a specialist, but perhaps this can help you.(It is a 
kind of "Bunch-Kaufmann-Partlett-algorithm)..
Cheers,
Stef.
=============================================================
#include <octave/oct.h>
 
#include <iostream.h>
 
#include <octave/pager.h>
#include <octave/symtab.h>
#include <octave/variables.h>
 
#include <octave/f77-fcn.h>
#include <octave/lo-error.h>
 
//#include <octave/mx-inlines.cc>
#include "BKPsolve_fine.h"
 
extern "C" {
 int F77_FCN (dsysvx, DSYSVX) (const char&,const char&,const int&,const 
int&,const double*,
                               const int&,double*,const int&,int*,const 
double*,const int&,
                               double*,const int&,double&,double*,double*,
                               double*,const int&,int*,int&);
 }
DEFUN_DLD (BKPsolve_fine, args, nargout,
  "[X,INFO]=BKPsolve_fine(A,B)\n\
\n\
Geeft X, waar A*X=B, opgelost met het\n\
Bunch-Kauffman-Partlett-algorithme voor \n\
niet positief definiete, symmetrische stelsels.
(iteratief verfijnd)
(.oct)")
{
  int nargin = args.length();
  if (nargin != 2 || nargout > 2)
    {print_usage ("BKPsolve_fine");
     return octave_value();
    }
  Matrix A = args(0).matrix_value();
  ColumnVector B = args(1).column_vector_value();
  int nrA=A.rows();
  int ncA=A.columns();
  int nrB=B.length();
  int arg_is_empty = empty_arg("BKPsolve",nrA,ncA);
  if (arg_is_empty<0)
    return octave_value();
  else if (arg_is_empty>0)
    return octave_value_list(3,Matrix());
  if (nrA!=ncA)
    {gripe_square_matrix_required ("BKPsolve");
     return octave_value();
    }
  if (nrA!=nrB)
    {(*current_liboctave_error_handler) ("BKPsolve: A en B moeten evenveel 
rijen tellen");
     return octave_value ();
    }
  char fact='N';
  char uplo='U';
  int n=nrA;
  int nrhs=1;
  double *pA = A.fortran_vec();
  int lda=n;
  int ldaf=n;
  Array2<double> af(ldaf,n);
  double *paf=af.fortran_vec();
  Array<int> ipiv(n);
  int *pipiv=ipiv.fortran_vec();
  double *pB=B.fortran_vec();
  int ldb=n;
 int ldx=n;
  Matrix X(ldx,nrhs);
  double *pX=X.fortran_vec();
  double rcond;
  Array<double> ferr(nrhs);
  double *pferr=ferr.fortran_vec();
  Array<double> berr(nrhs);
  double *pberr=berr.fortran_vec();
  int lwork=64*n;
  ColumnVector work(lwork);
  double *pwork=work.fortran_vec();
  Array<int> iwork(n);
  int *piwork=iwork.fortran_vec();
  int info=0;
 
 
  F77_XFCN 
(dsysvx,DSYSVX,(fact,uplo,n,nrhs,pA,lda,paf,ldaf,pipiv,pB,ldb,pX,ldx,
                          rcond,pferr,pberr,pwork,lwork,piwork,info));
 
  cout << "RCOND:       " << rcond << endl;
  if (info != 0)
   {
    cout << "LET OP: DSYSVX_error: INFO = " << info << endl;
                          rcond,pferr,pberr,pwork,lwork,piwork,info));
 
  cout << "RCOND:       " << rcond << endl;
  if (info != 0)
   {
    cout << "LET OP: DSYSVX_error: INFO = " << info << endl;
   }
  else
   {
/*    int lwork_opt=int(work(0));
    if (lwork != lwork_opt)
     {
      cout<<"LWORK        = "<<lwork<<endl;
      cout<<"LWORK optimal= "<<lwork_opt<<endl;
     }
*/
   }
  octave_value_list retval;
  if (nargout ==  1) retval(0)=X;
  else
    {retval(1)=double(info);retval(0)=X;}
  return retval;
}



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