He solucionado el problema anterior, era un problema al hacer el printf de la matriz, y ahora me calcula toda la matriz transpuesta, con las columnas y las filas bien. Pero el contenido de las celdas no es bueno. Empieza rellenando los valores de la fila1 con los valores de la columna1 de la matriz original, pero en el 4º o el 5º valor cambia de columna y empieza a coger los valores de la columna 2 de la matriz original... No se si me he explicado bien :s
Agradeceria mucho si alguien me pudiera ayudar, ya que es para un proyecto de clase. Este el codigo corregido:
Código c:
Ver original#include <stdio.h>
#include <QtCore/QCoreApplication>
#include <iostream>
#include <cmath>
using namespace std;
void affichage_matrice(double**matrice, int n);
double* algo_solution (double **U, double**L, double*b, int n);
double** cholesky (double **matrice, int taille);
double* fcourbe(double *x,double *y,double *v,int longueur);
double* gnewton(double *x,double *y, double *v, int longueur, int longueur1);
double** transposee(double **chol,int col, int fil);
int main(int argc, char *argv[]){
QCoreApplication a(argc, argv);
int longueur;
int longueur1;
double *x;
double *y;
double *v;
//Allocation de memoire dynamique
x
=(double*)malloc(sizeof(double)*longueur
); y
=(double*)malloc(sizeof(double)*longueur
); v
=(double*)malloc(sizeof(double)*longueur1
);
longueur=10;
longueur1=4;
x[0]=0.1;
x[1]=0.3;
x[2]=0.7;
x[3]=1.2;
x[4]=1.6;
x[5]=2.2;
x[6]=2.7;
x[7]=3.1;
x[8]=3.5;
x[9]=3.9;
y[0]=0.558;
y[1]=0.569;
y[2]=0.176;
y[3]=-0.207;
y[4]=-0.133;
y[5]=0.132;
y[6]=0.055;
y[7]=-0.09;
y[8]=-0.069;
y[9]=0.027;
v[0]=1;
v[1]=2;
v[2]=2;
v[3]=1;
double *model;
model=gnewton(x,y,v,longueur,longueur1);
int z;
printf("Le vecteur v optimise est:\n"); for(z=0;z<longueur1;z++){
}
return (int)a.exec;
}
double** cholesky (double **matrice, int taille){
int i, j, m;
double b, a, b1, a1;
double **chol;
chol
=(double**)malloc(sizeof(double*)*taille
);
for(i=0;i<taille;i++){
chol
[i
]=(double*)calloc(taille
,sizeof(double));
}
int k;
double *matest;
double *Test;
double Tf;
int cont=1;
int cont1=1;
matest
=(double*)malloc(sizeof(double)*taille
); Test
=(double*)malloc(sizeof(double)*taille
);
//Initialisation des tableaux
int z=0;
for(z=0;z<taille;z++){
matest[z]=1;
}
//On teste si la matrice entrée est symétrique définie positive
for (i=0;i<taille;i++){
for(j=0;j<taille;j++)
{
if (abs(matrice
[i
][j
]-matrice
[j
][i
])>0.01 ) cont=0;
}
}
if (!cont)
printf ("\n La matrice n'est pas symetrique!");
for(int i=0;i<taille;i++){
Test[i]=0;
for(int j=0;j<taille;j++){
Test[i]=Test[i]+(matest[i]*matrice[i][j]);
Tf=Tf+(Test[i]*matest[i]);
if(Tf<0)
cont1=0;
}
}
if (!cont1)
printf ("\n La matrice n'est pas definie positive!\n");
if(cont && cont1){
//On calcule les valeurs des cases de la matrice decomposée L
for(i=0;i<taille;i++){
b=0;
for(j=0;j<=i-1;j++){
b1=chol[i][j]*chol[i][j];
b=b+b1;
}
chol
[i
][i
]=sqrt(matrice
[i
][i
]-b
);
for(k=i+1;k<taille;k++){
a=0;
for(m=0;m<=k-1;m++){
a1=chol[i][m]*chol[k][m];
a=a+a1;
}
chol[k][i]=(matrice[k][i]-a)/chol[i][i];
}
}
}
return chol;
}
//ALGORITMO DE CALCULO DE LA TRANSPUESTA
double** transposee(double **chol,int col, int fil){
double **cholin;
cholin
= (double**)malloc(fil
* sizeof(*cholin
)); for(int m = 0; m < fil; m++)
{
cholin
[m
] = (double*)malloc(col
* sizeof(**cholin
)); }
/* Initialisation */
for (int i=0; i < col; i++){
for (int j=i; j < fil; j++){
cholin[i][j] = chol[j][i];
}
}
return cholin;
}
void affichage_matrice(double**matrice, int n)
{
for (int i=0;i<n;i++)
{
for (int j=0;j<n;j++)
{
cout<<matrice[i][j]<<" ";
}
cout<<endl;
}
cout<<endl;
}
double* algo_remontee (double**U, double*b, int n)
{
double*x;
x=new double[n];
double somme=0;
x[n-1]=b[n-1]/U[n-1][n-1];
for(int i=n-2;i>=0;i--)
{
somme=0;
for(int j=i+1;j<n;j++)
{
somme+=U[i][j]*x[j];
}
x[i]=(b[i]-somme)/U[i][i];
}
return x;
}
double* algo_descente (double**L, double*b, int n)
{
double*x;
x=new double[n];
double somme=0;
x[0]=b[0]/L[0][0];
for(int i=1;i<n;i++)
{
somme=0;
for(int j=0;j<i;j++)
{
somme+=L[i][j]*x[j];
}
x[i]=(b[i]-somme)/L[i][i];
}
return x;
}
double* algo_solution(double **U, double**L, double*b, int n)
{
// on a A, on décompose en LLt, L=L et Lt=U
double*y=new double[n];
y=algo_descente(L,b,n);
return algo_remontee(U,y,n);
}
double* fcourbe(double *x,double *y,double *v,int longueur){
double *Fc=new double[longueur];
for(int i=0; i<longueur;i++){
Fc
[i
]=y
[i
]-(v
[0]*exp(-v
[1]*x
[i
])*sin((v
[2]*x
[i
])+v
[3]));
}
return Fc;
}
double** jacobien_courbe(double *x,double *v,int longueur, int longueur1){
double **J;
J
= (double**)malloc(longueur
* sizeof(*J
)); for(int m = 0; m < longueur; m++)
{
J
[m
] = (double*)malloc(longueur1
* sizeof(**J
)); }
for(int i=0; i<longueur; i++){
J
[i
][0]= -exp(-v
[1]*x
[i
])*sin((v
[2]*x
[i
])+v
[3]); J
[i
][1]= x
[i
]*v
[0]*exp(-v
[1]*x
[i
])*sin((v
[2]*x
[i
])+v
[3]); J
[i
][2]= -x
[i
]*v
[0]*exp(-v
[1]*x
[i
])*cos((v
[2]*x
[i
])+v
[3]); J
[i
][3]= -v
[0]*exp(-v
[1]*x
[i
])*cos((v
[2]*x
[i
])+v
[3]); }
return J;
int e; int d;
for(e=0;e<longueur;e++){
for(d=0;d<longueur1;d++){
}
}
}
double* gnewton(double *x,double *y, double *v, int longueur, int longueur1){
double *F;
double **J;
double **A;
A
= (double**)malloc(longueur1
* sizeof(*A
)); for(int m = 0; m < longueur1; m++){
A
[m
] = (double*)malloc(longueur1
* sizeof(**A
)); }
double *B=new double[longueur1];
double *C=new double[longueur1];
double **L;
double **Lt;
double *S;
double **Jt;
int p;
do{
F=fcourbe(x,y,v,longueur);
J=jacobien_courbe(x,v,longueur,longueur1);
int e; int d;
printf("La matrice J est:\n");
for(e=0;e<longueur;e++){
for(d=0;d<longueur1;d++){
}
}
Jt=transposee(J,longueur1,longueur);
printf("La matrice Jt est:\n");
for(e=0;e<longueur1;e++){
for(d=0;d<longueur;d++){
}
}
for(int i=0;i<longueur1;i++){
for(int j=0;j<longueur1;j++){
A[i][j]=0;
for(p=0;p<longueur;p++){
A[i][j]=A[i][j]+(Jt[i][p]*J[p][j]);
}
}
}
for(int i=0;i<longueur1;i++){
B[i]=0;
for(int j=0;j<longueur;j++){
B[i]=B[i]+(Jt[i][j]*(F[j]));
}
C[i]=-B[i];
}
for(e=0;e<longueur1;e++){
for(d=0;d<longueur1;d++){
}
}
L=cholesky(A,longueur1);
Lt=transposee(L,longueur1,longueur1);
S=algo_solution(L,Lt,C,longueur1);
for(e=0;e<longueur1;e++){
}
for(int k=0; k<longueur1;k++){
v[k]=v[k]+S[k];
}
for(d=0;d<longueur1;d++){
}
}while(sqrt(B
[0]*B
[0]+B
[1]*B
[1]+B
[2]*B
[2]+B
[3]*B
[3])>(pow(10,-6)));
return v;
}
