LAPACKを用いたBスプライン
連立方程式部分
http://nalab.mind.meiji.ac.jp/~mk/labo/studying-C/intro-lapack-2004Nov/node3.html
http://www.cazoo.jp/calcsoft/library/lapack/clapack/usage/bandstrage/index.html
http://www016.upp.so-net.ne.jp/sk-peo/skpeo_n_bmatrix.html
これがわかりよい
http://www.cazoo.jp/calcsoft/library/lapack/clapack/usage/dgbsv/index.html
「dgbsv は、倍精度実数型データで与えられた (Double precision)、 一般の (正値性を持たない General)、帯行列 (Band matrix) を係数行列とする 連立1次方程式を解くサブルーチンである。」
下のルーチンは一般行列（GE）を使っている dgesvを使ってる。
#include<stdio.h> #include<stdlib.h> double baseN(int r,int i,double t,double *all_knots){ double w1=0.0,w2=0.0; // printf("r:%d all_knots:%lf all_knotsp:%lf t:%lf\n",r,all_knots[i],all_knots[i+1],t); if(r==1){ if(t >=all_knots[i] && t<=all_knots[i+1]) return 1.0; else return 0.0; } else{ if((all_knots[i+r]-all_knots[i+1])!=0){ w1=((all_knots[i+r]-t)/(all_knots[i+r]-all_knots[i+1])) * baseN(r-1,i+1,t,all_knots); } if((all_knots[i+r-1]-all_knots[i]) !=0 ){ w2=((t-all_knots[i])/(all_knots[i+r-1]-all_knots[i])) * baseN(r-1,i,t,all_knots); } // printf("r:%d all_knots:%lf %lf\n",r,all_knots[i],all_knots[i+r],w1+w2); return (w1+w2); } } double internal_baseN(int r,int i,double t,double *all_knots){ return baseN(r,i,t,all_knots); } double spline(double x,double *d,int r,double *all_knots,int n){ int j; double sum=0; for(j=0;j<n;j++){ sum=sum+d[j]*internal_baseN(r,j,x,all_knots); // printf("d_j%d %lf %lf\n",j,d[j],internal_baseN(r,j,x,all_knots)); } return sum; } int main(){ int i,j,k,jj; int iis,iif,jjs,jjf; int l; double* x; double* y; double** a; double dx; double* A; double* knots; double* all_knots; double* s_knots; double* f_knots; double* d; double q; double sum; int order=3; int rank=order+1; //int number_of_data=1000; int number_of_data=2000; int number_of_knots=550; //int number_of_coefficients=number_of_knots-2; int n_all_knots=number_of_knots+2*rank-2; long n=number_of_knots+rank-2; long* ipiv; long nrhs=1,lda=n,ldb=n,info; A=(double*)malloc(sizeof(double)*n*n); ipiv=(long*)malloc(sizeof(long)*n); x=(double*)malloc(sizeof(double)*number_of_data); y=(double*)malloc(sizeof(double)*number_of_data); knots=(double*)malloc(sizeof(double)*number_of_knots); d=(double*)malloc(sizeof(double)*n); all_knots=(double*)malloc(sizeof(double)*n_all_knots); s_knots=(double*)malloc(sizeof(double)*(rank-1)); f_knots=(double*)malloc(sizeof(double)*(rank-1)); //Initilaize for(i=0;i<number_of_data;i++){ x[i]=i; y[i]=i*i; } for(i=0;i<n*n;i++){ A[i]=0; } dx=(x[number_of_data-1]+0.0)/(number_of_knots-1); for(i=0;i<number_of_knots;i++){ knots[i]=dx*i; printf("%d %lf\n",i,knots[i]); } for(i=0;i<rank-1;i++){ s_knots[i]=knots[0]; } for(i=0;i<rank-1;i++){ f_knots[i]=knots[number_of_knots-1]; printf("fknots %lf\n",f_knots[i]); } for(i=0;i<rank-1;i++){ all_knots[i]=s_knots[i]+0.00000; } for(i=0;i<number_of_knots;i++){ all_knots[rank+i-1]=knots[i]; } for(i=0;i<rank-1;i++){ all_knots[rank-1+number_of_knots+i]=f_knots[i]; } for(i=0;i<n_all_knots;i++){ printf("allknots %d %lf\n",i,all_knots[i]); } for(i=0;i<n;i++){ d[i]=0; } jj=0; for(i=0;i<number_of_data;i++){ while(all_knots[jj+rank]<x[i]){ jj=jj+1; // printf("%d\n",j); } j=jj-1; // printf("i: j:%d %d %lf %lf \n",i,j,all_knots[j+rank],internal_baseN(4,j+4,x[i],all_knots)); // for(j=0;j<n;j++){ d[j]=d[j]+y[i]*internal_baseN(rank,j,x[i],all_knots); // } d[j+1]=d[j+1]+y[i]*internal_baseN(rank,j+1,x[i],all_knots); d[j+2]=d[j+2]+y[i]*internal_baseN(rank,j+2,x[i],all_knots); d[j+3]=d[j+3]+y[i]*internal_baseN(rank,j+3,x[i],all_knots); d[j+4]=d[j+4]+y[i]*internal_baseN(rank,j+4,x[i],all_knots); } for(i=0;i<n;i++){ for(j=0;j<n;j++){ sum=0; for(k=0;k<number_of_data;k++){ // if(1==1){ //ii= //jj= iis=j-rank; jjs=j-rank; iif=j+1; jjf=j+1; if(iis<0){ iis=0; } if(jjs<0){ jjs=0; } if(iif>number_of_knots-1){ iif=number_of_knots-1; } if(jjf>number_of_knots-1){ jjf=number_of_knots-1; } if(x[k]>=knots[iis] && x[k]<=knots[iif] && x[k]>=knots[jjs] && x[k]<=knots[jjf]){ sum=sum+internal_baseN(rank,i,x[k],all_knots)*internal_baseN(rank,j,x[k],all_knots); } } // a[i][j]=sum; A[i+n*j]=sum; } // printf("%d\n",i); } /* for(i=0;i<11;i++){ printf("i%d %lf\n",i,baseN(rank,i,all_knots[12],all_knots)); } dgesv_(&n,&nrhs,A,&lda,ipiv,d,&ldb,&info); for(i=0;i<n;i++){ printf("d%d %lf\n",i,d[i]); } */ dgesv_(&n,&nrhs,A,&lda,ipiv,d,&ldb,&info); for(i=0;i<number_of_data;i++){ printf("%lf %lf %lf %lf\n",x[i],spline(x[i],d,rank,all_knots,n),y[i],internal_baseN(rank,9,x[i],all_knots)); } }
gcc -O3 spline_given_knots.c -llapack
cとフォートランの配列の違いの注意
http://www.vi.cs.is.nagoya-u.ac.jp/watanabe-lab/graduates/01year/hiroyuki/lab/clapack/index.html
