indice di condizionamento

[valeria_tedesco]

Ciao ragazzi! Percaso sapete se esiste una routine di Lapack che calcola l'indice di condizionamento di una matrice?
Volevo anche porre un'altra domanda:
per caso come si comporta la routine dgesv quando la matrice dei coefficienti è mal-condizionata?
Grazie!

[buttari]

Valeria,
devi usare la routine dgesvx. Questa routine e' equivalente a dgesv ma ha delle feature avanzate:

1) ti restituisce una stima dell'indice di condizionamento
2) fa una equilibratura della matrice attraverso scaling
3) se la soluzione non e' soddisfacente fa iterative refinement per migliorarla
4) restituisce anche una stima dell'error bound

immagino che questo risponda a entrambe le tue domande.

PS
nota che nessuno degli sviluppatori di LAPACK parla italiano (a parte me che ormai sono un ex). Quindi se vuoi avere risposte piu' dettagliate ti consiglio di scrivere in inglese.

[valeria_tedesco]

Grazie mille! Proverò a fare così!

[valeria_tedesco]

Hi! I used dgesvx to solve a linear system.
For some inputs the system is solved saccesfully. But for others it give me this error:

** On entry to DLATRS parameter number 5 had an illegal value

what do mean it?

thanks.

[Julien Langou]

can you give the calling sequence you are using?
For which value of EQUED, TRANS, N, NRHS, etc. do you see this error message?
-j

[valeria_tedesco]

Hello! When I call this function sometimes the system is solved succesfully,
but with some inputs of n and m the program give me this error:

** On entry to DLATRS parameter number 5 had an illegal value

This is my function:

void approssimante(int n, int m, double *p, double *q, double *a)
{

// I need that the matrix will be equilibrated if necessary, then copied to af and factored.
int N=n+m; //N=n+m number of linear equations
double *b;
double *noti; //I must solve bx=noti

char fact[1];
fact[0] = 'E';
char trans[1];
trans[0] = 'N';
char equed[1];
equed[0] = 'R';
int nrhs = 1;

int info;
int lda = N;
int ldaf = N;
int ldb = N;
int ldx = N;
int iwork = 4 * N;
double rcond;


//allocations
double *af=(double *)malloc(sizeof(double)*N*N);
double *R=(double*)malloc(sizeof(double)*N);
double *c=(double*)malloc(sizeof(double)*N);
double *x=(double*)malloc(sizeof(double)*N);
int *ipiv=(int*)malloc(sizeof(int)*N);
double *work=(double*)malloc(sizeof(double)*4*N);
double *ferr=(double*)malloc(sizeof(double)*N);
double *berr=(double*)malloc(sizeof(double)*N);

//initializations
for (int i = 0; i < 4*N; i++){
work[i] = 0;
}

for (int i = 0; i < N; i++){
ipiv[i] = 0;
R[i] = 0;
c[i] = 0;
x[i] = 0;
ferr[i] = 0;
berr[i] = 0;
}

b=costruisci(n,m,a); //b is a matrix NxN, N=n+m
noti=costruisciNoti(a,N);

dgesvx(fact, trans, &N, &nrhs, b, &lda, af, &ldaf, ipiv, equed, R, c, noti, &ldb, x, &ldx, &rcond, ferr, berr, work, &iwork, &info);

if (info<0)
{
cout << "L'elemento "<< -info << " è sbagliato nella chiamata";
p=NULL;q=NULL;
}
else
if (info>0)
{
cout << "La matrice è singolare\n";
p=NULL;q=NULL;
}
else //sistema compatibile
{
p[0]=a[0];
for (i=1;i<=n;++i) //copy of the firsts n-1 elements of the solution in the array p
p[i]=x[i-1];

q[0]=1;
for (i=1;i<=m;++i) //copy of the restant elements in the array q
q[i]=x[i+n-1];
}

}

Thanks!valy.

[buttari]

Valeria,
ferr and berr should be of size nrhs and not of size n.
See if this fixed your problem.

Alfredo