[Top][All Lists]

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

Re: Fortran

From: Stef Pillaert
Subject: Re: Fortran
Date: Tue, 15 May 2001 17:03:26 +0200

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)..
#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/>
#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,
Geeft X, waar A*X=B, opgelost met het\n\
Bunch-Kauffman-Partlett-algorithme voor \n\
niet positief definiete, symmetrische stelsels.
(iteratief verfijnd)
  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;
  cout << "RCOND:       " << rcond << endl;
  if (info != 0)
    cout << "LET OP: DSYSVX_error: INFO = " << info << endl;
  cout << "RCOND:       " << rcond << endl;
  if (info != 0)
    cout << "LET OP: DSYSVX_error: INFO = " << info << endl;
/*    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;
  return retval;

Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:
How to fund new projects:
Subscription information:

reply via email to

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