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