/*===========================*/
/* gccSEEP4nod.c             */
/*===========================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#define NDmax 3000
#define NEmax 5000
#define MAmax 10
#define NTmax NDmax*1

void DNMAT_S4(int ne,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1],
           double dnx[5],double dny[5],double *detJ,double a,double b);
void IPOINT4(int kk,double *a,double *b);
void STIFF_S4(int ne,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1],
           double Akx[NEmax+1],double Aky[NEmax+1],double eak[5][5]);
void CALV_S4(int ne,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1],double hvec[NDmax+1],
          double Akx[NEmax+1],double Aky[NEmax+1],double vx[NEmax+1],double vy[NEmax+1]);
void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][5],double sm[5][5],double tsm[NTmax+1][NTmax+1]);
void CHBAND10(int n1,int m1,double array[NTmax+1][NTmax+1]);
void CHBAND11(int n1,int m1,double array[NTmax+1][NTmax+1],double vector[NTmax+1]);
int IBAND(int mm,double a[NTmax+1][NTmax+1],double fact0);
void BAND(int mm,int ib,double a[NTmax+1][NTmax+1]);

/*------------------------------------------------------------------------------*/
int main(int argc, char* argv[])
{
    char fnameR[50];
    char fnameW[50];

    strcpy(fnameR,argv[1]);
    strcpy(fnameW,argv[2]);

    main_part(fnameR,fnameW);

    return 0;
}
/*------------------------------------------------------------------------------*/

int main_part(char fnameR[50],char fnameW[50])
{
    /*２次元飽和定常浸透流解析*/
    FILE *fin,*fout;
    char dat[256]="";
    char sv[20]="";
    int i,j,k,n,ia,ja,ne;
    char strcom[256]="";    /* 書出用コメント */
    int idan;               /* 0：水平断面，1：鉛直断面 */
    int NODT;               /* 節点総数 */
    int NELT;               /* 要素総数 */
    int MATEL;              /* 材料種類数 */
    int KOH;                /* 全水頭指定節点数 */
    int KOQ;                /* 流量指定節点数 */
    int nt;                 /* 全自由度（総節点数×１[１節点自由度]）*/
    int mm;                 /* 温度固定点処理後自由度 */
    int ib;                 /* バンド幅 */
    double elarea;          /* 要素半面積 */
    double fact0=1.0E-30;   /* 0判定 */
    clock_t t_sta,t_end;
    double cal_time;
    time_t jikoku;
    struct tm *lt;
    char *p;
    /*-----------------------------------------------------*/
    /*配列寸法宣言*/
    /*-----------------------------------------------------*/
    int    kakom[NEmax+1][5];    /*要素構成接点番号*/
    double Akx[NEmax+1];         /*要素x方向透水係数*/
    double Aky[NEmax+1];         /*要素y方向透水係数*/
    int    matno[NEmax+1];       /*材料種別No*/
    double wm[NEmax+1][3];       /*作業用材料物性記憶*/
    double x[NDmax+1];           /*節点x座標*/
    double y[NDmax+1];           /*節点y座標*/
    int    nokh[NDmax+1];        /*全水頭指定節点番号*/
    int    nokq[NDmax+1];        /*流量指定節点番号*/
    double Hinp[NDmax+1];        /*指定全水頭入力*/
    double Qinp[NDmax+1];        /*指定流量入力*/
    double vx[NEmax+1];          /*x方向流速ベクトル*/
    double vy[NEmax+1];          /*y方向流速ベクトル*/
    int    index[NTmax+1];       /*温度固定点処理用*/
    static double tak[NTmax+1][NTmax+1];/*全体透水行列*/
    double hvec[NTmax+1];        /*全体全水頭ベクトル*/
    double qvec[NTmax+1];        /*全体節点流量ベクトル*/
    double reac[NTmax+1];        /*作業用ベクトル*/
    static double wak[NTmax+1][NTmax+1];/*水頭指定節点対応列記憶*/
    double eak[5][5];            /*要素透水行列*/

    /*時間計測開始*/
    t_sta=clock();
    /******************************************************/
    /*データ入力*/
    /******************************************************/
    fin=fopen(fnameR,"r");

printf("%s\n",fnameR);
    /*--------------------------------------------------------------------------------------*/
    /*書出用コメント入力*/
    /*--------------------------------------------------------------------------------------*/
    fgets(dat,sizeof dat,fin);strncpy(strcom,dat,strlen(dat)-1);
    /*--------------------------------------------------------------------------------------*/
    /*解析断面(0：水平，1：鉛直)，節点数，要素数，材料種類数,全水頭固定節点数，流量指定節点数*/
    /*--------------------------------------------------------------------------------------*/
    fgets(dat,sizeof dat,fin);
    idan  =atoi(strtok(dat,","));
    NODT  =atoi(strtok(NULL,","));
    NELT  =atoi(strtok(NULL,","));
    MATEL =atoi(strtok(NULL,","));
    KOH   =atoi(strtok(NULL,","));
    KOQ   =atoi(strtok(NULL,","));
    nt=NODT*1; /* 全自由度 */

printf("%d %d %d %d %d %d\n",idan,NODT,NELT,MATEL,KOH,KOQ);

    /*--------------------------------------------------------------------------------------*/
    /*材料特性（Akx,Aky）*/
    /*--------------------------------------------------------------------------------------*/
    for(i=1;i<=MATEL;i++){
        fgets(dat,sizeof dat,fin);
        wm[i][1]=atof(strtok(dat,","));
        wm[i][2]=atof(strtok(NULL,","));
    }
    /*--------------------------------------------------------------------------------------*/
    /*要素構成節点と材料種別（node-1,node-2,node-3,matno）*/
    /*--------------------------------------------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        fgets(dat,sizeof dat,fin);
        kakom[ne][1]=atoi(strtok(dat,","));
        kakom[ne][2]=atoi(strtok(NULL,","));
        kakom[ne][3]=atoi(strtok(NULL,","));
        kakom[ne][4]=atoi(strtok(NULL,","));
        matno[ne]=atoi(strtok(NULL,","));
        for(i=1;i<=MATEL;i++){
            if(i==matno[ne]){
                Akx[ne]=wm[i][1];
                Aky[ne]=wm[i][2];
            }
        }
    }
    /*--------------------------------------------------------------------------------------*/
    /*節点座標（x,y）*/
    /*--------------------------------------------------------------------------------------*/
    for(i=1;i<=NODT;i++){
        fgets(dat,sizeof dat,fin);
        x[i]     =atof(strtok(dat,","));
        y[i]     =atof(strtok(NULL,","));
    }
    /*--------------------------------------------------------------------------------------*/
    /*全水頭指定節点番号，指定全水頭*/
    /*--------------------------------------------------------------------------------------*/
    if(0<KOH){
        for(i=1;i<=KOH;i++){
            fgets(dat,sizeof dat,fin);
            nokh[i]=atoi(strtok(dat,","));
            Hinp[i]=atof(strtok(NULL,","));
        }
    }
    /*--------------------------------------------------------------------------------------*/
    /*流量指定節点番号，指定流量*/
    /*--------------------------------------------------------------------------------------*/
    if(0<KOQ){
        for(i=1;i<=KOQ;i++){
            fgets(dat,sizeof dat,fin);
            nokq[i]=atoi(strtok(dat,","));
            Qinp[i]=atof(strtok(NULL,","));
        }
    }
    fclose(fin);
    /******************************************************/
    /* 入力確認出力 */
    /******************************************************/
    fout=fopen(fnameW,"w");
    fprintf(fout,"%s\n",strcom);
    fprintf(fout,"idan,NODT,NELT,MATEL,KOH,KOQ\n");
    sprintf(sv,"%d",idan)   ;strcpy(dat,sv);
    sprintf(sv,"%d",NODT)   ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",NELT)   ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",MATEL)  ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOH)    ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOQ)    ;strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
    fprintf(fout,"*input-data\n");
    fprintf(fout,"node,x,y,Hinp,Qinp\n");
    for(i=1;i<=NODT;i++){
        sprintf(sv,"%d",i)           ;strcpy(dat,sv);
        sprintf(sv,"%e",x[i])        ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",y[i])        ;strcat(dat,",");strcat(dat,sv);
        strcpy(sv,"--");
        for(j=1;j<=KOH;j++){if(i==nokh[j])sprintf(sv,"%e",Hinp[j]);}
                                      strcat(dat,",");strcat(dat,sv);
        strcpy(sv,"--");
        for(j=1;j<=KOQ;j++){if(i==nokq[j])sprintf(sv,"%e",Qinp[j]);}
                                      strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fprintf(fout,"element,node-1,node-2,node-3,Akx,Aky,matno\n");
    for(ne=1;ne<=NELT;ne++){
        sprintf(sv,"%d",ne)           ;strcpy(dat,sv);
        sprintf(sv,"%d",kakom[ne][1]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",kakom[ne][2]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",kakom[ne][3]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",kakom[ne][4]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",Akx[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",Aky[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",matno[ne])    ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fclose(fout);
    /*要素面積確認と異常終了通知*/
    for(ne=1;ne<=NELT;ne++){
        i=kakom[ne][1];j=kakom[ne][2];k=kakom[ne][3];
        elarea=0.5*((x[k]-x[j])*y[i]+(x[i]-x[k])*y[j]+(x[j]-x[i])*y[k]);
        if(elarea<fact0){
            /* 結果出力（Appendモード）*/
            fout=fopen(fnameW,"a+");
            fprintf(fout,"area of element is less or equal to 0\n");
            fclose(fout);
            printf("要素面積<=0：異常終了通知");
            fclose(fout);
            return 0;
        }
        i=kakom[ne][3];j=kakom[ne][4];k=kakom[ne][1];
        elarea=0.5*((x[k]-x[j])*y[i]+(x[i]-x[k])*y[j]+(x[j]-x[i])*y[k]);
        if(elarea<fact0){
            /* 結果出力（Appendモード）*/
            fout=fopen(fnameW,"a+");
            fprintf(fout,"area of element is less or equal to 0\n");
            fclose(fout);
            printf("要素面積<=0：異常終了通知");
            fclose(fout);
            return 0;
        }
    }
    /*******************************************************/
    /*メイン処理計算*/
    /*******************************************************/
    /*-----------------------------------------------------*/
    /*全体行列クリア*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){
        reac[i]=0.0;
        hvec[i]=0.0;
        qvec[i]=0.0;
        for(j=1;j<=nt;j++){
            tak[i][j]=0.0;
    }}
    /*-----------------------------------------------------*/
    /*全体行列組立*/
    /*-----------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        STIFF_S4(ne,kakom,x,y,Akx,Aky,eak);
        ASMAT(ne,4,1,kakom,eak,tak);
    }
    /*-----------------------------------------------------*/
    /*全水頭指定節点対応列記憶，前処理*/
    /*-----------------------------------------------------*/
    if(0<KOH){
        for(k=1;k<=KOH;k++){
            for(i=1;i<=nt;i++){
                wak[i][k]=tak[i][nokh[k]];
        }}
        for(k=1;k<=KOH;k++){
            for(i=1;i<=nt;i++){
                hvec[i]=hvec[i]+Hinp[k]*wak[i][k];
        }}
    }
    if(0<KOQ){
        for(k=1;k<=KOQ;k++){
            qvec[nokq[k]]=Qinp[k];
        }
    }
    for(i=1;i<=nt;i++){
        qvec[i]=qvec[i]-hvec[i];
    }
    /*-----------------------------------------------------*/
    /*全水頭指定節点処理*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){index[i]=i;}
    if(0<KOH){
        for(i=1;i<=KOH;i++){
            n=nokh[i];
            index[n]=0;
        }
    }
    mm=0;
    for(i=1;i<=nt;i++){
        if(index[i]>0){
            mm=mm+1;
            index[mm]=i;
        }
    }
    for(i=1;i<=mm;i++){
        ia=index[i];
        reac[i]=qvec[ia];
        for(j=1;j<=mm;j++){
            ja=index[j];
            tak[i][j]=tak[ia][ja];
    }}
    /*-----------------------------------------------------*/
    /*バンド幅計算とフルマトリックスのバンド化記憶*/
    /*-----------------------------------------------------*/
    ib=IBAND(mm,tak,fact0);
    BAND(mm,ib,tak);
    /*-----------------------------------------------------*/
    /*対角項確認と異常終了通知*/
    /*-----------------------------------------------------*/
    for(i=1;i<=mm;i++){
        if(fabs(tak[i][1])<fact0){
            fout=fopen(fnameW,"a+");
            fprintf(fout,"diagonal value is less or equal to 0\n");
            fclose(fout);
            printf("対角項<=0：異常終了通知");
            fclose(fout);
            return 0;
        }
    }
    /*-----------------------------------------------------*/
    /*コレスキー法*/
    /*-----------------------------------------------------*/
    CHBAND10(mm,ib,tak);/*三角行列*/
    CHBAND11(mm,ib,tak,reac);
    /*-----------------------------------------------------*/
    /*全水頭算定*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){hvec[i]=0.0;}
    for(i=1;i<=mm;i++){ia=index[i];hvec[ia]=reac[i];}
    for(k=1;k<=KOH;k++){hvec[nokh[k]]=Hinp[k];}
    /*-----------------------------------------------------*/
    /*節点流量算定*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){qvec[i]=0.0;}
    if(0<KOQ){
        for(k=1;k<=KOQ;k++){
            qvec[nokq[k]]=Qinp[k];
        }
    }
    if(0<KOH){
        for(k=1;k<=KOH;k++){
            for(i=1;i<=nt;i++){
                qvec[nokh[k]]=qvec[nokh[k]]+wak[i][k]*hvec[i];
        }}
    }
    /*-----------------------------------------------------*/
    /*要素平均流速算定*/
    /*-----------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        CALV_S4(ne,kakom,x,y,hvec,Akx,Aky,vx,vy);
    }
    /*******************************************************/
    /*結果出力*/
    /*******************************************************/
    /*結果出力（Appendモード）*/
    fout=fopen(fnameW,"a+");
    /*節点値*/
    fprintf(fout,"*output-data\n");
    fprintf(fout,"node,x,y,hvec,qvec,P-head\n");
    /*節点番号,x座標,y座標,全水頭,節点流量,圧力水頭*/
    /*節点流量は，流入が正・流出が負の値*/
    for(i=1;i<=NODT;i++){
        sprintf(sv,"%d",i)       ;strcpy(dat,sv);
        sprintf(sv,"%e",x[i])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",y[i])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",hvec[i]) ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",qvec[i]) ;strcat(dat,",");strcat(dat,sv);
        strcpy(sv,"--");
        if(idan==1)sprintf(sv,"%e",hvec[i]-y[i]);
                                  strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fprintf(fout,"element,vx,vy\n");
    /*流速ベクトル*/
    for(ne=1;ne<=NELT;ne++){
        sprintf(sv,"%d",ne)     ;strcpy(dat,sv);
        sprintf(sv,"%e",vx[ne])  ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",vy[ne])  ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    /*計算完了*/
    t_end=clock();
    cal_time=(double)(t_end-t_sta)/CLOCKS_PER_SEC;
    time(&jikoku);
    lt=localtime(&jikoku);
    p=asctime(lt);
    fprintf(fout,"NODT=%d nt=%d mm=%d ib=%d\n",NODT,nt,mm,ib);
    fprintf(fout,"Calculation time=%.3f(sec)\n",cal_time);
    fprintf(fout,"%s\n",p);
    fclose(fout);

    printf("NODT=%d nt=%d mm=%d ib=%d\n",NODT,nt,mm,ib);
    printf("計算完了\n");
    printf("計算時間：%.3f sec\n",cal_time);
    printf("%s\n",p);
    return 0;

}
/*------------------------------------------------------------------------------*/

void DNMAT_S4(int ne,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1],
           double dnx[5],double dny[5],double *detJ,double a,double b)
{
    int i,j,k,l;
    double J11,J12,J21,J22;
    double dn1a,dn2a,dn3a,dn4a;
    double dn1b,dn2b,dn3b,dn4b;

    i=kakom[ne][1];
    j=kakom[ne][2];
    k=kakom[ne][3];
    l=kakom[ne][4];

    /*[dN/dadN/db]*/
    dn1a=-0.25*(1.0-b);
    dn2a= 0.25*(1.0-b);
    dn3a= 0.25*(1.0+b);
    dn4a=-0.25*(1.0+b);
    dn1b=-0.25*(1.0-a);
    dn2b=-0.25*(1.0+a);
    dn3b= 0.25*(1.0+a);
    dn4b= 0.25*(1.0-a);
    /*Jacobiマトリックスとdet(J)*/
    J11=dn1a*x[i]+dn2a*x[j]+dn3a*x[k]+dn4a*x[l];
    J12=dn1a*y[i]+dn2a*y[j]+dn3a*y[k]+dn4a*y[l];
    J21=dn1b*x[i]+dn2b*x[j]+dn3b*x[k]+dn4b*x[l];
    J22=dn1b*y[i]+dn2b*y[j]+dn3b*y[k]+dn4b*y[l];
    *detJ=J11*J22-J12*J21;
    /*[dN/dxdN/dy]*/
    dnx[1]=( J22*dn1a-J12*dn1b)/ *detJ;
    dnx[2]=( J22*dn2a-J12*dn2b)/ *detJ;
    dnx[3]=( J22*dn3a-J12*dn3b)/ *detJ;
    dnx[4]=( J22*dn4a-J12*dn4b)/ *detJ;
    dny[1]=(-J21*dn1a+J11*dn1b)/ *detJ;
    dny[2]=(-J21*dn2a+J11*dn2b)/ *detJ;
    dny[3]=(-J21*dn3a+J11*dn3b)/ *detJ;
    dny[4]=(-J21*dn4a+J11*dn4b)/ *detJ;
}
/*------------------------------------------------------------------------------*/

void IPOINT4(int kk,double *a,double *b)
{
    /*面積積分におけるGauss積分点*/
    *a=0.0;
    *b=0.0;
    switch(kk){
        case 1: *a=-0.5773502692; *b=-0.5773502692;break;
        case 2: *a= 0.5773502692; *b=-0.5773502692;break;
        case 3: *a= 0.5773502692; *b= 0.5773502692;break;
        case 4: *a=-0.5773502692; *b= 0.5773502692;break;
    }
}
/*------------------------------------------------------------------------------*/

void STIFF_S4(int ne,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1],
           double Akx[NEmax+1],double Aky[NEmax+1],double eak[5][5])
{
    /*要素透水行列*/
    int kk,i,j;
    double dnx[5];
    double dny[5];
    double detJ,a,b;
    double wx[5][5];
    double wy[5][5];

    for(i=1;i<=4;i++){
        for(j=1;j<=4;j++){
            eak[i][j]=0.0;
    }}

    for(kk=1;kk<=4;kk++){
        IPOINT4(kk,&a,&b);
        DNMAT_S4(ne,kakom,x,y,dnx,dny,&detJ,a,b);
        for(i=1;i<=4;i++){
            for(j=1;j<=4;j++){
                wx[i][j]=0.0;
                wy[i][j]=0.0;
        }}
        for(i=1;i<=4;i++){
            for(j=1;j<=4;j++){
                wx[i][j]=wx[i][j]+dnx[i]*dnx[j];
                wy[i][j]=wy[i][j]+dny[i]*dny[j];
        }}
        for(i=1;i<=4;i++){
            for(j=1;j<=4;j++){
                eak[i][j]=eak[i][j]+(Akx[ne]*wx[i][j]+Aky[ne]*wy[i][j])*detJ;
        }}
    }
}
/*------------------------------------------------------------------------------*/

void CALV_S4(int ne,int kakom[NEmax+1][5],double x[NDmax+1],double y[NDmax+1],double hvec[NDmax+1],
          double Akx[NEmax+1],double Aky[NEmax+1],double vx[NEmax+1],double vy[NEmax+1])
{
    /*流速ベクトル算定*/
    int kk,j;
    int inod[5];
    double dnx[5];
    double dny[5];
    double detJ,a,b,wvx,wvy;

    vx[ne]=0.0;
    vy[ne]=0.0;
    for(j=1;j<=4;j++){
        inod[j]=kakom[ne][j];
    }
    for(kk=1;kk<=4;kk++){
        IPOINT4(kk,&a,&b);
        DNMAT_S4(ne,kakom,x,y,dnx,dny,&detJ,a,b);
        wvx=0.0;
        wvy=0.0;
        for(j=1;j<=4;j++){
            wvx=wvx+dnx[j]*hvec[inod[j]];
            wvy=wvy+dny[j]*hvec[inod[j]];
        }
        vx[ne]=vx[ne]+wvx;
        vy[ne]=vy[ne]+wvy;
    }
    vx[ne]=-Akx[ne]*vx[ne]*0.25;
    vy[ne]=-Aky[ne]*vy[ne]*0.25;
}
/*------------------------------------------------------------------------------*/

void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][5],double sm[5][5],double tsm[NTmax+1][NTmax+1])
{
    /*剛性行列組立[nod：1要素節点数，nhen：1節点の自由度]*/
    int i,j,k,l,ki,kj,ks,js,kik,kjl,ksk,jsl;

    for(i=1;i<=nod;i++){
        for(j=1;j<=nod;j++){
            ki=(kakom[ne][i]-1)*nhen;
            kj=(kakom[ne][j]-1)*nhen;
            ks=(i-1)*nhen;
            js=(j-1)*nhen;
            for(k=1;k<=nhen;k++){
                for(l=1;l<=nhen;l++){
                    kik=ki+k;
                    kjl=kj+l;
                    ksk=ks+k;
                    jsl=js+l;
                    tsm[kik][kjl]=tsm[kik][kjl]+sm[ksk][jsl];
                }
            }
        }
    }

}
/*------------------------------------------------------------------------------*/

void CHBAND10(int n1,int m1,double array[NTmax+1][NTmax+1])
{
    /*三角行列*/
    int i,j,k,mm,mn;
    double z;

    array[1][1]=sqrt(array[1][1]);
    for(j=2;j<=m1;j++){
        array[1][j]=array[1][j]/array[1][1];
    }
    for(i=2;i<=n1;i++){
        if(n1-i+1>m1){
            mm=m1;
        }else{
            mm=n1-i+1;
        }
        for(j=1;j<=mm;j++){
            z=array[i][j];
            if(j!=m1){
                if(m1-j+1<i){
                    mn=m1-j+1;
                }else{
                    mn=i;
                }
                for(k=2;k<=mn;k++){
                    z=z-array[i-k+1][k]*array[i-k+1][j+k-1];
                }
                if(j>1){
                    array[i][j]=z/array[i][1];
                }else{
                    array[i][1]=sqrt(z);
                }
            }else{
                array[i][j]=z/array[i][1];
            }
        }
    }
}
/*--------------------------------------------------------------------------------*/

void CHBAND11(int n1,int m1,double array[NTmax+1][NTmax+1],double vector[NTmax+1])
{
    /*前進消去・後退代入*/
    int i,j,mm;
    double z;

    vector[1]=vector[1]/array[1][1];
    for(i=2;i<=n1;i++){
        z=vector[i];
        if(i>m1){
            mm=m1;
        }else{
            mm=i;
        }
        for(j=2;j<=mm;j++){
            z=z-array[i-j+1][j]*vector[i-j+1];
        }
        vector[i]=z/array[i][1];
    }

    vector[n1]=vector[n1]/array[n1][1];
    for(i=n1-1;i>=1;i--){
        z=vector[i];
        for(j=2;j<=m1;j++){
            if(i+j-1>n1)break;
            z=z-array[i][j]*vector[i+j-1];
        }
        vector[i]=z/array[i][1];
    }
}
/*--------------------------------------------------------------------------------*/

int IBAND(int mm,double array[NTmax+1][NTmax+1],double fact0)
{
    int ib,i,j;
    /*-----------------------------------------------------*/
    /*バンド幅計算
    /*-----------------------------------------------------*/
    ib=1;
    for(i=1;i<=mm-1;i++){
        for(j=mm;j>=i+1;j--){
            if(fact0<fabs(array[i][j]))break;
        }
        if(ib<=j-i+1)ib=j-i+1;
    }
    if(mm<ib)ib=mm;
    return ib;
}
/*--------------------------------------------------------------------------------*/

void BAND(int mm,int ib,double array[NTmax+1][NTmax+1])
{
    int i,j,k,mn;
    double work[mm+1];
    /*-----------------------------------------------------*/
    /*フルマトリックスのバンド化記憶
    /*-----------------------------------------------------*/
    for(i=2;i<=mm;i++){
        if(i+ib-1<mm){
            mn=ib;
        }else{
            mn=mm-i+1;
        }
        for(j=1;j<=ib;j++){
            work[j]=0.0;
        }
        for(k=1;k<=mn;k++){
            work[k]=array[i][i+k-1];
        }
        for(j=1;j<=ib;j++){
            array[i][j]=work[j];
        }
    }
}
/*--------------------------------------------------------------------------------*/