[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
octave+arpack=problem
From: |
ctkaczyk-oct |
Subject: |
octave+arpack=problem |
Date: |
Thu, 17 Jun 2004 12:21:24 +0200 |
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
-------------------------------------------------------------
- octave+arpack=problem,
ctkaczyk-oct <=