[Top][All Lists]

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

Re: octave+arpack=problem

From: Paul Kienzle
Subject: Re: octave+arpack=problem
Date: Thu, 17 Jun 2004 07:11:41 -0400

Your Fortran calls need character length arguments.  See
F77_*CHAR* details on the following wiki page:

Also, rather than using p = malloc(n*m*sizeof(double)), you
can use:

        Matrix P(m,n);
        double* p = P.fortran_vec();

This makes sure you are using whatever memory scheme
Octave is using, and also cleans up memory for you when
you are done.

Even better, use:


which also cleans up after itself.

Paul Kienzle

On Jun 17, 2004, at 6:21 AM, address@hidden wrote:


I've post this to address@hidden
and was asked to post it here, so:


I am trying to implement function eigs(), as Matlab does it, using an Arpack library. Unfortunately I?ve stacked at the beginning and can't move any step

Problem is as follows...
Simple program using arpack written in C works OK. But almost the same code used as a function in oct-file produces wrong results. For example aprack
function asked for starting vector always returns [NaN NaN NaN ....]
Of course, the same problem arise in C++ instead if C.

I think of some memory problem. It looks like aprack functions are referring to wrong place in memory when called from octave function. However, some of them do it all right. I've debug Arpack some how, and notice that i.e. this starting
vector looks all right till it is normalized.

I wonder if:
1)anybody has faced similar problem?
2)anybody is interested in helping me? If so, I'll put some more details.

Czarek Tkaczyk

And here comes details:
I've put DEFUN_DLD below, but I haven't written meanings of parameters used in dsaupd(). I suppose it doesn't matter, because almost the same code compiled with gcc, instead of mkoctfile, works OK.

#include <oct.h>
#include "f77-fcn.h"

#define I(i)      ((i)-1)

extern "C"
       int F77_FUNC (dsaupd, DSAUPD)(int& ido, char* bmat, int& n,
                                     char* which, int& nev, double& tol,
                                     double* resid, int& ncv, double* v,
                                     int& ldv, int* iparam, int* ipntr,
                                     double* workd, double* workl, int& lworkl,
                                     int& info );

DEFUN_DLD(dssimp, args, , "the simplest possible use of arpack function dsaupd()")
  octave_value_list retval;

  int maxn = 256, maxnev = 10, maxncv = 25;
  int ldv = maxn;

//    Code written in this manner:
//      Matrix v (ldv , maxncv);
//      fortran_function ( ..., v.fortran_vec (), ...);
//    generates the same problem

   double* v, *workl, *workd, *d, *resid, *ax;
   v     = (double*)malloc(ldv * maxncv *sizeof(double));
   workl = (double*)malloc(maxncv * (maxncv + 8) *sizeof(double));
   workd = (double*)malloc(3*maxn*sizeof(double));
   d     = (double*)malloc(maxncv * 2*sizeof(double));
   resid = (double*)malloc(maxn*sizeof(double));
   ax    = (double*)malloc(maxn*sizeof(double));

   bool   select[maxncv];
   int    iparam[11],
          ipntr [11];

   int    ido, n, nev, ncv, lworkl, info, ierr,
          j, nx, ishfts, maxitr, mode1, nconv;
   bool   rvec;
   double tol, sigma;

   nx = 7;
   n  = nx * nx;

   nev   = 4;
   ncv   = 20;
   char bmat[2]  = "I";
   char which[3] = "LM";

   lworkl = ncv*(ncv+8);
   tol = 0.1;
   info = 0;
   ido = 0;

   ishfts = 1;
   maxitr = 300;
   mode1 = 1;
   iparam[I(1)] = ishfts;
   iparam[I(3)] = maxitr;
   iparam[I(7)] = mode1;

//   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//    And here come problems.
//    dsaupd() works like this:
//      * generate starting vector
//      * do the first step of Lanczos factorisation
// After debugging I've notice, that starting vector looks like this:
//     [NaN NaN NaN ...]
//    In this case factorisation returns error
//   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        F77_XFCN (dsaupd, DSAUPD, ( ido, bmat, n, which, nev, tol, resid,
                                    ncv, v, ldv, iparam, ipntr, workd, workl,
                                    lworkl, info ) );
//    Now, resid = [NaN NaN NaN ...]
        if (f77_exception_encountered)
           error ("unrecoverable error in DSAUPD");
           return retval;

  return retval;


 Czarek Tkaczyk

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

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

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]