【C言語】LAPACKを使った多変数の最小二乗法の計算方法【Cygwin】

Cygwinを用いて、多変数の最小二乗法を行ったので、備忘録として、下記に留める。

となるを求めるプログラム。（今回は5次まで）
の順でスペースで区切られた近似したいデータファイル(test.dat)を読み込んで計算します。
gcc sample.c -llapack -lblas -lm のコマンドでコンパイルします。

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

int dim,dim1,dim2,dim3;
void matrix_transposition(int,int,double  **A,double  **t_A);//行列転置
void matrix_initialize(int,int,double  **A); //行列要素すべて０
void matrix_inverse(int,double **A,double **B);
void matrix_inverse2(int,double **A,double **B);
void matrix_calculate(int,int,int, double **A,double **B,double **C);
double g_calc(int,double T,double omega);
void test_printf(int,int,double **A,char *);//テスト用表示

#define SIZE 21
#define NSIZE 30

int main(void){
  int i,j;
  FILE *fp;
  int ret;
  int cnt = 0;
  
 
  char fname[100] = "test.dat";
	//printf("filename=");
   //scanf("%s",&fname);
	
	
  fp = fopen( fname , "r" );
  if( fp == NULL ){
   printf( "%sファイルが開けません\n", fname );
    return -1;
  }

  while((ret = getc(fp)) != EOF) {
		 if(ret == '\n') cnt++;
  }
 
	
  //int NSIZE = cnt;
	printf("NSIZE = %d\n",NSIZE);
 static double S_total,S[NSIZE],SI[NSIZE],sigma,rS[SIZE];
	
	double **G;
	G = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		G[i] = malloc(sizeof(double) * NSIZE);
	}
    double **t_G;
	t_G = malloc(sizeof(double *) * NSIZE);
	for (i=0;i<NSIZE;i++) {
		t_G[i] = malloc(sizeof(double) * SIZE);
	}
	double **temp;
	temp = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		temp[i] = malloc(sizeof(double) * SIZE);
	}
	double **i_temp;
	i_temp = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		i_temp[i] = malloc(sizeof(double) * SIZE);
	}
	double **temp2;
	temp2 = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		temp2[i] = malloc(sizeof(double) * 1);
	}
	double **I;
	I = malloc(sizeof(double *) * NSIZE);
	for (i=0;i<NSIZE;i++) {
		I[i] = malloc(sizeof(double) * 1);
	}
	double *T;
	T = malloc(sizeof(double ) * NSIZE);
	double *omega;
	omega = malloc(sizeof(double ) * NSIZE);
	
	double **a;
	a = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		a[i] = malloc(sizeof(double) * 1);
	}
	double **unit;
	unit = malloc(sizeof(double *) * SIZE);
	for (i=0;i<SIZE;i++) {
		unit[i] = malloc(sizeof(double) * SIZE);
	}
  rewind( fp );
	
  for(i=0;i<NSIZE;i++){
    ret = fscanf( fp, "%lf %lf %lf", &I[i][0],&T[i],&omega[i]) ;
  	//T[i] = T[i];
  	//omega[i] = omega[i];
  	//printf("%e %e %e\n",I[i][0],T[i],omega[i]);
  }
  fclose(fp);
	
	for(i=0;i<SIZE;i++){
		for(j=0;j<NSIZE;j++){
			G[i][j] = g_calc(i,T[j],omega[j]);
			//printf("G[%d][%d] =%e\n",i,j,G[i][j]);
		}
	}
	//test_printf(SIZE,NSIZE,G,"G");
	matrix_transposition(SIZE,NSIZE,G,t_G);
	//test_printf(SIZE,NSIZE,t_G,"t_G");
	matrix_calculate(SIZE,NSIZE,SIZE,G,t_G,temp);
	test_printf(SIZE,SIZE,temp,"temp");
	matrix_inverse(SIZE,temp,i_temp);
	//matrix_calculate(SIZE,SIZE,SIZE,i_temp,temp,unit);
	//test_printf(SIZE,SIZE,unit,"unit");
	//puts("pass3");
	matrix_calculate(SIZE,NSIZE,1,G,I,temp2);
	//test_printf(SIZE,SIZE,i_temp,"i_temp");
	//test_printf(SIZE,SIZE,i_temp,"i_temp");
	matrix_calculate(SIZE,SIZE,1,i_temp,temp2,a);
	
	for(i=0;i<SIZE;i++){
		printf("a[%d] = %e\n",i,a[i][0]);
	}
	
	
	for(j=0;j<NSIZE;j++){
		SI[j]=0.;
		for(i=0;i<SIZE;i++){
			SI[j] += a[i][0]* g_calc(i,T[j],omega[j]);
		}
		S[j] = (I[j][0] - SI[j])*(I[j][0] - SI[j])/((double) NSIZE);
		rS[j] = sqrt(S[j]*((double) NSIZE));
		printf("%e\n",rS[j]);
	}
	
	S_total = 0.;
	for(j=0;j<NSIZE;j++){
		S_total += S[j]; 
	}
	
	
	sigma = sqrt(S_total);
	printf("standard deviation = %e",sigma);
	
	
	free( G );	
	free( t_G );	
	free( temp );	
	free( i_temp );	
	free( temp2 );	
	free( I );	
	free( T );	
	free( omega );
	free( a );	
}

double g_calc(int i,double T,double omega){
	double ans;
	switch (i) {
	case 0:
		ans =1.0;
		return ans;
		break;

	case 1:
		ans = T;
		return ans;
		break;
		
	case 2:
		ans = omega;
		return ans;
		break;
		
	case 3:
		ans =T*T;
		return ans;
		break;
		
	case 4:
		ans =T*omega;
		return ans;
		break;
		
	case 5:
		ans =omega*omega;
		return ans;
		break;
		
	case 6:
		ans = T*T*T;
		return ans;
		break;

	case 7:
		ans = T*T*omega;
		return ans;
		break;
		
	case 8:
		ans = T*omega*omega;
		return ans;
		break;
		
	case 9:
		ans =omega*omega*omega;
		return ans;
		break;
		
		
		
	case 10:
		ans =T*T*T*T;
		return ans;
		break;
		
	case 11:
		ans = T*T*T*omega;
		return ans;
		break;

	case 12:
		ans = T*T*omega*omega;
		return ans;
		break;
		
	case 13:
		ans = T*omega*omega*omega;
		return ans;
		break;
		
	case 14:
		ans =omega*omega*omega*omega;
		return ans;
		break;
		
		
	case 15:
		ans =T*T*T*T*T;
		return ans;
		break;
		
	case 16:
		ans = T*T*T*T*omega;
		return ans;
		break;

	case 17:
		ans = T*T*T*omega*omega;
		return ans;
		break;
		
	case 18:
		ans = T*T*omega*omega*omega;
		return ans;
		break;
		
	case 19:
		ans =T*omega*omega*omega*omega;
		return ans;
		break;
		
	case 20:
		ans =omega*omega*omega*omega*omega;
		return ans;
		break;	
		
	default:
		break;
	}
}


void matrix_calculate(int dim1,int dim2,int dim3, double **A,double **B,double **C){
	int i,j,k;
	matrix_initialize(dim1,dim3,C);
	for(i=0;i<dim1;i++){
		for(j=0;j<dim3;j++){
	  		for(k=0;k<dim2;k++){
				C[i][j]=C[i][j]+A[i][k]*B[k][j];
	  		}
	 	}
	}
}


void matrix_inverse2(int dim,double **A,double **B){
	
	double buf; //一時的なデータを蓄える
	int i,j,k; //カウンタ
	int n=dim;  //配列の次数
	double inv_a[dim][dim];
 	double C[dim][dim];
	
	for(i=0;i<n;i++){
		for(j=0;j<n;j++){
 			C[i][j]= A[i][j];
 		}
	}
	
	//単位行列を作る
	for(i=0;i<n;i++){
 		for(j=0;j<n;j++){
 			if(i==j){ B[i][j]=1.0;}
 			else{B[i][j] = 0.0;}
 		}
	}
	//掃き出し法
	for(i=0;i<n;i++){
		buf=1/C[i][i];
		for(j=0;j<n;j++){
			C[i][j]*=buf;
 			B[i][j]*=buf;
		}
		for(j=0;j<n;j++){
 			if(i!=j){
			buf=C[j][i];
			for(k=0;k<n;k++){
				C[j][k]-=C[i][k]*buf;
				B[j][k]-=B[i][k]*buf;
  			}
			}
		}
	}
}

void matrix_inverse(int dim,double **A,double **B){
  	int i,j;
  	long   m       = dim ; // 行のサイズ
  	long   n       = dim ; // 列のサイズ
  	long   lda     = dim ; // mと同じ値
	double *C;
	C = malloc(sizeof(double) * dim);
	
 	for(i=0;i<dim;i++){
	 for(j=0;j<dim;j++){
		C[dim*i+j] = A[i][j];
	 }
	}
	//test_printf(dim,dim,B,"B");
  	long   info           ; 
  	long   ipiv[dim]     ;
  	long   lwork   = dim ;
  	double work[dim]     ;

  	dgetrf_( &m, &n, C, &lda, ipiv, &info);
	if(info !=0){puts("error_LU");}
  	dgetri_( &n, C, &lda, ipiv, work, &lwork, &info);
	if(info !=0){puts("error_inverse");}

	for(i=0;i<dim;i++){
		for(j=0;j<dim;j++){
			B[i][j] = C[dim*i+j];
		}
	}
}

//行列要素すべて0
void matrix_initialize(int dim1,int dim2,double **A){
	int i,j;
	for(i=0;i<dim1;i++){
		for(j=0;j<dim2;j++){
			A[i][j]=0.0;
		}
	}
}

//行列転置
void matrix_transposition(int dim1,int dim2,double **A,double **t_A){
	int i,j;
	for(i=0;i<dim1;i++){
		for(j=0;j<dim2;j++){
			t_A[j][i]=A[i][j];
		}
	}
}

//テスト用表示
void test_printf(int dim1,int dim2,double **A,char *name){
	int i,j;
	for(i=0;i<dim1;i++){
		for(j=0;j<dim2;j++){
			printf("%s[%d,%d] = %e \n",name,i,j,A[i][j]);
		}	
	}
	printf("\n");
}