[C] Interpolazione mediante spline cubiche

[salvodel]

Salve a tutti dovrei scrivere un funzione che mi interpoli una serie di punti con delle spline cubiche. Al momento sto partendo da zero e volevo chiedervi alcune cose:
1 - dove posso trovare in rete la parte matematica relativa alle spline?:wtf:
2 - siccome da quel che ho capito dovrò invertire una matrice, come posso fare? Ci sono delle librerie che facciano questa operazione in modo ottimizzato?
Grazie a tutti

[wizard1993]

2 - siccome da quel che ho capito dovrò invertire una matrice, come posso fare? Ci sono delle librerie che facciano questa operazione in modo ottimizzato?
Grazie a tutti

le intel math kernel ma costano (per win)
http://www.intel.com/cd/software/products/asmo-na/eng/307757.htm
per lin sono free

[salvodel]

le intel math kernel ma costano (per win)
http://www.intel.com/cd/software/products/asmo-na/eng/307757.htm
per lin sono free

Grazie mille per la risposta ma ahi me avevo dato gia un occhiata a questo thread
http://www.hwupgrade.it/forum/showthread.php?t=1751833.

Cercando ho trovato qualcosa per C# ma quello che mi stupisce è la necessita di dover scrivere cose basilari. Mi spiego: io adesso ho dovuto implementare una funzione per interpolare dei punti con una spline del 3 ordine e va bene ma al suo interno dovendo invertire una matrice cosa faccio?
1 - Inverto la matrice in base alle regole dell'algebra matriciale?
2 - Implemento io il metodo di Gauss-Jordan?
Possibile che non ci siano nella rete queste operazioni semplici gia scritte? La soluzione uno è banale ma all'aumentare della complessità della matrice perdo parecchio tempo di calco.
Grazie

[salvodel]

Ho trovato una funzione che fa l'inversione con il metodo di Gauss-Jordan solo che non mi funziona.
Posto il listato che ho trovato con la parte che ho inserito io, se qualcuno riesce a buttargli un occhio giusto per capire dove sbaglio. In pratica come risultato non mi da X ma B.:muro:
Grazie

/******************************************************************************/
/* Perform Gauss-Jordan elimination with row-pivoting to obtain the solution to
* the system of linear equations
* A X = B
*
* Arguments:
* lhs - left-hand side of the equation, matrix A
* rhs - right-hand side of the equation, matrix B
* nrows - number of rows in the arrays lhs and rhs
* ncolsrhs- number of columns in the array rhs
*
* The function uses Gauss-Jordan elimination with pivoting. The solution X to
* the linear system winds up stored in the array rhs; create a copy to pass to
* the function if you wish to retain the original RHS array.
*
* Passing the identity matrix as the rhs argument results in the inverse of
* matrix A, if it exists.
*
* No library or header dependencies, but requires the function swaprows, which
* is included here.
*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// swaprows - exchanges the contents of row0 and row1 in a 2d array
void swaprows(double** arr, long row0, long row1) {
double* temp;
temp=arr[row0];
arr[row0]=arr[row1];
arr[row1]=temp;
}

// gjelim
void gjelim(double** lhs, double** rhs, long nrows, long ncolsrhs) {

// augment lhs array with rhs array and store in arr2
double** arr2=new double*[nrows];
for (long row=0; row<nrows; ++row)
arr2[row]=new double[nrows+ncolsrhs];

for (long row=0; row<nrows; ++row) {
for (long col=0; col<nrows; ++col) {
arr2[row][col]=lhs[row][col];
}
for (long col=nrows; col<nrows+ncolsrhs; ++col) {
arr2[row][col]=rhs[row][col-nrows];
}
}

// perform forward elimination to get arr2 in row-echelon form
for (long dindex=0; dindex<nrows; ++dindex) {
// run along diagonal, swapping rows to move zeros in working position
// (along the diagonal) downwards
if ( (dindex==(nrows-1)) && (arr2[dindex][dindex]==0)) {
return; // no solution
} else if (arr2[dindex][dindex]==0) {
swaprows(arr2, dindex, dindex+1);
}
// divide working row by value of working position to get a 1 on the
// diagonal
if (arr2[dindex][dindex] == 0.0) {
return;
} else {
double tempval=arr2[dindex][dindex];
for (long col=0; col<nrows+ncolsrhs; ++col) {
arr2[dindex][col]/=tempval;
}
}

// eliminate value below working position by subtracting a multiple of
// the current row
for (long row=dindex+1; row<nrows; ++row) {
double wval=arr2[row][dindex];
for (long col=0; col<nrows+ncolsrhs; ++col) {
arr2[row][col]-=wval*arr2[dindex][col];
}
}
}

// backward substitution steps
for (long dindex=nrows-1; dindex>=0; --dindex) {
// eliminate value above working position by subtracting a multiple of
// the current row
for (long row=dindex-1; row>=0; --row) {
double wval=arr2[row][dindex];
for (long col=0; col<nrows+ncolsrhs; ++col) {
arr2[row][col]-=wval*arr2[dindex][col];
}
}
}

// assign result to replace rhs
for (long row=0; row<nrows; ++row) {
for (long col=0; col<ncolsrhs; ++col) {
rhs[row][col]=arr2[row][col+nrows];
}
}

for (long row=0; row<nrows; ++row)
delete[] arr2[row];
delete[] arr2;
}

int main()
{
int i,j;
long dim;
double **A,**b;

dim=3;

A = (double **)malloc(dim*sizeof(double *));
if(A==NULL)
printf("MEmoria esaurita!\n");
for(i=0; i<dim; i++)
{
A[i]=(double *)malloc(dim*sizeof(double));
if(A==NULL)
printf("Memoria esaurita!\n");
}

b = (double **)malloc(dim*sizeof(double *));
if(b==NULL)
printf("MEmoria esaurita!\n");
for(i=0; i<dim; i++)
{
b[i]=(double *)malloc(1*sizeof(double));
if(b==NULL)
printf("Memoria esaurita!\n");
}

for(i=0;i<dim;i++)
for(j=0;j<dim;j++)
{
printf("A[%d][%d]= ",i+1,j+1);
scanf("%lf",&A[i][j]);
}
for(i=0;i<dim;i++)
{
for(j=0;j<dim;j++)
printf("A[%d][%d]= %3.1lf\t",i+1,j+1,A[i][j]);
printf("\n");
}
for(i=0;i<dim;i++)
{
printf("b[%d][1]=",i+1);
scanf("%lf",&b[i][1]);
}
for(i=0;i<dim;i++)
printf("b[%d][1]= %3.1lf\n",i+1,b[i][1]);

gjelim(A, b, dim, 1);

printf("\n--------------\n\nSoluzione\n\n");
for(i=0;i<dim;i++)
printf("b[%d][1]= %3.1lf\n",i+1,b[i][1]);

system("PAUSE");

}
:help::ave:

[Netskate]

se vuoi io ho una implementazione di gauss jordan funzionante. però è fatta in fortran :D

[salvodel]

se vuoi io ho una implementazione di gauss jordan funzionante. però è fatta in fortran :D

Grazie mille per la proposta ma mi sa che faccio prima a riprendere il libro di calocolo numero e a riscrivermi l'algoritmo...non sono molto esperto di fortran.:(
Casomai parto da questo che ho trovato e cerco la magagna. Ciao

[Netskate]

come vuoi. puoi anche provare a cercare la libreria imsl che contiene un sacco di funzioni matematiche. fra cui anche il calcolo dell'inversa

[salvodel]

come vuoi. puoi anche provare a cercare la libreria imsl che contiene un sacco di funzioni matematiche. fra cui anche il calcolo dell'inversa

Grazie mille per il suggerimento:mano: .