octave+arpack=problem

[ctkaczyk-oct]

Hello,

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

> Hello,
>
> 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
> forward.
>
> 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.
>
> Regards,
> 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;

}




 Regards,
 Czarek Tkaczyk



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