gslで線形方程式、連立一次方程式をとく。

gslで線形方程式、連立一次方程式をとく。

#include<stdio.h>
#include<gsl/gsl_linalg.h>

void solve(double** m,double* v,int len){
        using namespace std;
        int s;
        double *mb;
        int k=0;
        mb=(double*)malloc(sizeof(double)*len*len);
        for(int i=0;i<len;i++){
                for(int j=0;j<len;j++){
                      
                        mb[k]=m[i][j];
                        k++;

                }

        }


        gsl_matrix_view m2=gsl_matrix_view_array(mb,len,len);
        gsl_vector_view v2=gsl_vector_view_array(v,len);
        gsl_vector *x=gsl_vector_alloc(len);
        gsl_permutation *p =gsl_permutation_alloc(len);
        gsl_linalg_LU_decomp(&m2.matrix,p,&s);
        gsl_linalg_LU_solve(&m2.matrix,p,&v2.vector,x);
        printf("x=\n");

        gsl_vector_fprintf(stdout,x,"%g");
        gsl_permutation_free(p);
        gsl_vector_free(x);

//      return *v;

}


int main(){

        double** m;
        double* v;
        double* x;
        int len=3;


        m=(double**) malloc(sizeof(double*)*(len));
        v=(double*) malloc(sizeof(double)*(len));
        x=(double*) malloc(sizeof(double)*(len));

        for(int i=0;i<len;i++){
                m[i]=(double*)malloc(sizeof(double)*(len));
        }

        m[0][0]=3;
        m[0][1]=2;
        m[0][2]=5;
        m[1][0]=7;
        m[1][1]=4;
        m[1][2]=1;
        m[2][0]=9;
        m[2][1]=2;
        m[2][2]=6;

        v[0]=2;
        v[1]=4;
        v[2]=7;

        solve(m,v,3);


}

x=
0.818182
マイナス0.454545
0.0909091

帰り値あり

#include<stdio.h>
#include<gsl/gsl_linalg.h>

void solve(double** m,double* v,int len,double* output){
        using namespace std;
        int s;
        double *mb;
        int k=0;
        mb=(double*)malloc(sizeof(double)*len*len);
        for(int i=0;i<len;i++){
                for(int j=0;j<len;j++){
//                      printf("%lf\n",m[i][j]);
                        mb[k]=m[i][j];
                        k++;

                }

        }


        gsl_matrix_view m2=gsl_matrix_view_array(mb,len,len);
        gsl_vector_view v2=gsl_vector_view_array(v,len);
        gsl_vector *x=gsl_vector_alloc(len);
        gsl_permutation *p =gsl_permutation_alloc(len);
        gsl_linalg_LU_decomp(&m2.matrix,p,&s);
        gsl_linalg_LU_solve(&m2.matrix,p,&v2.vector,x);
        printf("x=\n");

        gsl_vector_fprintf(stdout,x,"%g");
        gsl_permutation_free(p);
        //gsl_vector_free(x);
        for(int i=0;i<len;i++){
                output[i]=gsl_vector_get(x,i);

        };
//      return *v;

}


int main(){

        double** m;
        double* v;
        double* x;
        double* output;
        int time_length=3;


        m=(double**) malloc(sizeof(double*)*(time_length));
        v=(double*) malloc(sizeof(double)*(time_length));
        x=(double*) malloc(sizeof(double)*(time_length));

        output=(double*) malloc(sizeof(double)*(time_length));
        for(int i=0;i<time_length;i++){
                m[i]=(double*)malloc(sizeof(double)*(time_length));
        }

        m[0][0]=3;
        m[0][1]=2;
        m[0][2]=5;
        m[1][0]=7;
        m[1][1]=4;
        m[1][2]=1;
        m[2][0]=9;
        m[2][1]=2;
        m[2][2]=6;


/*
        m[0][0]=1;
        m[0][1]=0;
        m[0][2]=0;
        m[1][0]=0;
        m[1][1]=1;
        m[1][2]=0;
        m[2][0]=0;
        m[2][1]=0;
        m[2][2]=1;
*/
        v[0]=2;
        v[1]=4;
        v[2]=7;

        solve(m,v,3,output);
        printf("answer\n");
        for(int i;i<time_length;i++){
                printf("%lf\n",output[i]);

        }

}

出力

x=
0.818182

• 0.454545

0.0909091
answer
0.818182

• 0.454545

0.090909