[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?
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]

Originariamente inviato da salvodel

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/pro...eng/307757.htm
per lin sono free

[salvodel]

Originariamente inviato da wizard1993

le intel math kernel ma costano (per win)
http://www.intel.com/cd/software/pro...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.
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");

}

[Netskate]

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

[salvodel]

Originariamente inviato da Netskate

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

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]

Originariamente inviato da 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

Grazie mille per il suggerimento .