/*===========================*/
/* gccFRAME_GAL.c            */
/*===========================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#define NDmax 500       /*最大節点数*/
#define NEmax 500       /*最大要素数*/
#define MAmax 10        /*最大材料種別数*/
#define NTmax NDmax*3+1 /*最大自由度数*/
int main_part(char fnameR[50],char fnameW[50],int kprog,double coefAB,double deltaS0,double coefS,int maxload);
double CAL_LAM01(int nt,double deltaS,double Spara,double dis0[NTmax+1],double df[NTmax+1]);
double CAL_LAM02(int nt,double deltaS,double coefAB,double *Spara,double dis0[NTmax+1],double df[NTmax+1]);
double CAL_LAMNEW(int nt,double Spara,double dis0[NTmax+1],double disr[NTmax+1],double df[NTmax+1]);
void DATAOUT(char fnameW[50],int NODT,int NELT,
    int iload,int itaration,int nocon,double deltaS,double lam,double dtest,double ftest,
    double x[NDmax+1],double y[NDmax+1],double dist[NTmax+1],
    double ff[NTmax+1],double ftvec[NTmax+1],double reforce[NEmax+1][7]);
void FRAMESM(int ne,int kakom[NEmax+1][3],double x[NDmax+1],double y[NDmax+1],double dist[NTmax+1],
    double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double reforce[NEmax+1][7],double sm[7][7]);
void ELMSM(int ne,double xx1,double yy1,double xx2,double yy2,
    double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double sme[7][7]);
void CALST(int ne, int kakom[NEmax+1][3], double x[NDmax+1],double y[NDmax+1],double dis[NTmax+1],
    double dist[NTmax+1],double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double fevec[7],double reforce[NEmax+1][7]);
void ROTATE(int ne, int kakom[NEmax+1][3], double x[NDmax+1],double y[NDmax+1],
    double dis[NTmax+1],double dist[NTmax+1],double *deltai,double *deltaj);
void HENMX(int ne,double xx1,double yy1,double xx2,double yy2,double hen[7][7]);
void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][3],double sm[7][7],double tsm[NTmax+1][NTmax+1]);
void ASVEC(int ne,int nod,int nhen,int kakom[NEmax+1][3],double fevec[7],double ftvec[NTmax+1]);
int 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 array[NTmax+1][NTmax+1],double fact0);
void BAND(int mm,int ib,double array[NTmax+1][NTmax+1]);
int MATINV(int n,double ary[NTmax+1][NTmax+1]);
/*--------------------------------------------------------------------------------*/
int main(int argc, char* argv[])
{
    char fnameR[50]; /*入力ファイル名*/
    char fnameW[50]; /*出力ファイル名*/
    int kprog;       /*第一中間点設定方法(1:方法1，2：方法2)*/
    double coefAB;   /*スケーリングパラメータ算定用係数*/
    double deltaS0;  /*弧長増分法パラメータ(初期弧長)*/
    double coefS;    /*弧長に乗じる係数*/
    int maxload;     /*荷重増分数最大値*/

    strcpy(fnameR,argv[1]);
    strcpy(fnameW,argv[2]);
    kprog   =atoi(argv[3]);
    coefAB  =atof(argv[4]);
    deltaS0 =atof(argv[5]);
    coefS   =atof(argv[6]);
    maxload =atoi(argv[7]);

    main_part(fnameR,fnameW,kprog,coefAB,deltaS0,coefS,maxload);

/*
    printf("リターンキーを押してください");
    getchar();
*/
    return 0;
 
}
/*--------------------------------------------------------------------------------*/
int main_part(char fnameR[50],char fnameW[50],int kprog,double coefAB,double deltaS0,double coefS,int maxload)
{
    /*微小変形有限変位平面骨組構造解析*/
    FILE *fin,*fout;
    char dat[1000]="";
    char sv[50]="";
    int i,j,k,l,m,ne,n,ia,ja;
    char strcom[50]="";   /*コメント*/
    int NODT;             /*節点総数*/
    int NELT;             /*要素総数*/
    int MATEL;            /*材料種別数*/
    int KOX;              /*x方向変位既知節点数*/
    int KOY;              /*y方向変位既知節点数*/
    int KOZ;              /*回転変位既知節点数*/
    int NF;               /*載荷節点数*/
    int nt;               /*全自由度*/
    int mm;               /*変位既知節点処理語の自由度*/
    int ib;               /*バンド幅*/
    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][3];            /*要素構成節点番号*/
    double Em[NEmax+1];                  /*要素弾性係数*/
    double AA[NEmax+1];                  /*要素断面積*/
    double AI[NEmax+1];                  /*要素断面二次モーメント*/
    int    matno[NEmax+1];               /*材料種別No*/
    double wm[MAmax+1][4];               /*材料種別作業用*/
    double x[NDmax+1];                   /*初期節点x座標*/
    double y[NDmax+1];                   /*初期節点y座標*/
    int    nokx[NDmax+1];                /*x方向変位既知節点番号*/
    int    noky[NDmax+1];                /*y方向変位既知節点番号*/
    int    nokz[NDmax+1];                /*回転変位既知節点番号*/
    int    index[NTmax+1];               /*変位既知節点処理のためのインデックス*/
    static double tsm[NTmax+1][NTmax+1]; /*全体剛性マトリックス*/
    static double twk[NTmax+1][NTmax+1]; /*全体剛性マトリックス*/
    double df[NTmax+1];                  /*全体節点外力ベクトル基準増分(入力値)*/
    double ff[NTmax+1];                  /*全体節点外力ベクトル(荷重)*/
    double ftvec[NTmax+1];               /*全体節点力ベクトル*/
    double dis[NTmax+1];                 /*全体節点変位増分*/
    double dist[NTmax+1];                /*全体節点変位*/
    double dsav[NTmax+1];                /*全体節点変位記憶(釣合点)*/
    double dis0[NTmax+1];                /*全体節点変位増分(基準外力相当増分)*/
    double disr[NTmax+1];                /*全体節点変位増分(不平衡力相当増分)*/
    double reac[NTmax+1];                /*作業用ベクトル*/
    double work[NTmax+1];                /*作業用ベクトル*/
    double sm[7][7];                     /*要素剛性マトリックス*/
    double fevec[7];                     /*要素節点荷重ベクトル*/
    double hen[7][7];                    /*座標変換行列*/
    double dtest,dtest0;                 /*変位収束判定用変数*/
    double ftest,ftest0;                 /*不平衡力収束判定用変数*/
    int maxita=50;                       /*繰返数最大値*/
    int itaration;                       /*収束計算回数インクリメント*/
    int iload;                           /*荷重増分インクリメント*/
    int nocon;                           /*基準回数での収束か否かのフラグ*/
    double reforce[NEmax+1][7];          /*要素断面力*/
    double refsave[NEmax+1][7];          /*要素断面力記憶(釣合点)*/
    double lam,dlam;                     /*弧長増分法パラメータ*/
    double deltaS;                       /*弧長増分法パラメータ(計算に用いる弧長)*/
    double Spara;                        /*弧長増分法パラメータ(スケーリングパラメータ)*/
    double fjdg0,fjdg;                   /*弧長変更用*/
    int jdg,jdginv;                      /*特異行列判定*/
    double fend0,fend1;                  /*荷重の動きの監視*/
    double dend0,dend1;                  /*変位の動きの監視*/
    int jflagc;                          /*計算終了フラグ*/

    /*時間計測開始*/
    t_sta=clock();
    /*******************************************************/
    /*データ入力*/
    /*******************************************************/
    fin=fopen(fnameR,"r");
    /*--------------------------------------------------------------------------------------*/
    /*書出用コメント入力*/
    /*--------------------------------------------------------------------------------------*/
    fgets(dat,sizeof dat,fin);strncpy(strcom,dat,strlen(dat)-1);
    /*--------------------------------------------------------------------------------------*/
    /*節点数，要素数，材料種別数][X方向変位既知節点数，Y方向変位既知節点数，回転変位既知節点数，載荷点数*/
    /*--------------------------------------------------------------------------------------*/
    fgets(dat,sizeof dat,fin);
    NODT   =atoi(strtok(dat,","));
    NELT   =atoi(strtok(NULL,","));
    MATEL  =atoi(strtok(NULL,","));
    KOX    =atoi(strtok(NULL,","));
    KOY    =atoi(strtok(NULL,","));
    KOZ    =atoi(strtok(NULL,","));
    NF     =atoi(strtok(NULL,","));
    nt=NODT*3;
    /*------------------------------------------------------------------------------------------*/
    /*材料特性（Em,AA,AI）*/
    /*------------------------------------------------------------------------------------------*/
    for(i=1;i<=MATEL;i++){
        fgets(dat,sizeof dat,fin);
        wm[i][1]=atof(strtok(dat,","));
        wm[i][2]=atof(strtok(NULL,","));
        wm[i][3]=atof(strtok(NULL,","));
    }
    /*------------------------------------------------------------------------------------------*/
    /*要素構成節点と材料種別（node-1,node-2,matno）*/
    /*------------------------------------------------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        fgets(dat,sizeof dat,fin);
        kakom[ne][1]=atoi(strtok(dat,","));
        kakom[ne][2]=atoi(strtok(NULL,","));
        matno[ne]   =atoi(strtok(NULL,","));
        for(i=1;i<=MATEL;i++){
            if(i==matno[ne]){
                Em[ne]   =wm[i][1];
                AA[ne]   =wm[i][2];
                AI[ne]   =wm[i][3];
            }
        }
    }
    /*--------------------------------------------------------------------------------------*/
    /*節点座標（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<KOX){
        for(i=1;i<=KOX;i++){
            fgets(dat,sizeof dat,fin);
            nokx[i]=atoi(strtok(dat,","));
        }
    }
    if(0<KOY){
        for(i=1;i<=KOY;i++){
            fgets(dat,sizeof dat,fin);
            noky[i]=atoi(strtok(dat,","));
        }
    }
    if(0<KOZ){
        for(i=1;i<=KOZ;i++){
            fgets(dat,sizeof dat,fin);
            nokz[i]=atoi(strtok(dat,","));
        }
    }
    /*--------------------------------------------------------------------------------------*/
    /*節点番号,x方向荷重,y方向荷重,モーメント荷重*/
    /*--------------------------------------------------------------------------------------*/
    for(i=1;i<=nt;i++){df[i]=0.0;}
    if(0<NF){
        for(i=1;i<=NF;i++){
            fgets(dat,sizeof dat,fin);
            n=atoi(strtok(dat,","));
            df[3*n-2]=atof(strtok(NULL,","));
            df[3*n-1]=atof(strtok(NULL,","));
            df[3*n]  =atof(strtok(NULL,","));
        }
    }
    fclose(fin);
    /*******************************************************/
    /*入力確認出力*/
    /*******************************************************/
    fout=fopen(fnameW,"w");
    fprintf(fout,"%s\n",strcom);
    fprintf(fout,"NODT,NELT,MATEL,KOX,KOY,KOZ,NF\n");
    sprintf(sv,"%d",NODT)     ;strcpy(dat,sv);
    sprintf(sv,"%d",NELT)     ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",MATEL)    ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOX)      ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOY)      ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",KOZ)      ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%d",NF)       ;strcat(dat,",");strcat(dat,sv);
    sprintf(sv,"%.3f",deltaS0);strcat(dat,",");strcat(dat,sv);
    fprintf(fout,"%s\n",dat);
    fprintf(fout,"*input-data\n");
    fprintf(fout,"node,x,y,dfx,dfy,dfz,x-fix,y-fix,z-fix\n");
    for(i=1;i<=NODT;i++){
        k=0;l=0;m=0;
        if(0<KOX){for(j=1;j<=KOX;j++){if(i==nokx[j])k=1;}}
        if(0<KOY){for(j=1;j<=KOY;j++){if(i==noky[j])l=1;}}
        if(0<KOZ){for(j=1;j<=KOZ;j++){if(i==nokz[j])m=1;}}
        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",df[3*i-2])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",df[3*i-1])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",df[3*i])      ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",k)           ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",l)           ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",m)           ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fprintf(fout,"element,node-1,node-2,E,A,I,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,"%e",Em[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",AA[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",AI[ne])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",matno[ne])    ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fclose(fout);
    /*******************************************************/
    /*メイン処理計算*/
    /*******************************************************/
    /*-----------------------------------------------------*/
    /*変位既知節点自由度記憶*/
    /*-----------------------------------------------------*/
    for(i=1;i<=nt;i++){index[i]=i;}
    if(0<KOX){for(i=1;i<=KOX;i++){n=3*nokx[i]-2;index[n]=0;}}
    if(0<KOY){for(i=1;i<=KOY;i++){n=3*noky[i]-1;index[n]=0;}}
    if(0<KOZ){for(i=1;i<=KOZ;i++){n=3*nokz[i]  ;index[n]=0;}}
    mm=0;
    for(i=1;i<=nt;i++){if(index[i]>0){mm=mm+1;index[mm]=i;}}
    /*-----------------------------------------------------*/
    /*ベクトルクリア*/
    /*-----------------------------------------------------*/
    for(ne=1;ne<=NELT;ne++){
        for(i=1;i<=6;i++){
            reforce[ne][i]=0.0;
            refsave[ne][i]=0.0;
    }}
    for(i=1;i<=nt;i++){
        ff[i]=0.0;
        dist[i]=0.0;
        dsav[i]=0.0;
        dis0[i]=0.0;
        disr[i]=0.0;
        reac[i]=0.0;
    }
    /*-----------------------------------------------------*/
    /*方法1におけるスケーリングパラメータ設定*/
    /*-----------------------------------------------------*/
    ftest=0.0;
    for(i=1;i<=nt;i++){
        ftest=ftest+df[i]*df[i];
    }
    Spara=sqrt(coefAB/ftest);
    ftest=0.0;
    /*******************************************************/
    /*荷重増分過程*/
    /*******************************************************/
    jflagc=0;/*正常終了フラグ初期設定*/
    iload=0;
    deltaS=deltaS0;
    do{
        iload=iload+1;
        jdg=0;
again:
        itaration=0;
        nocon=0;
        lam=0.0;
        dtest0=0.0;
        ftest0=0.0;
        for(ne=1;ne<=NELT;ne++){for(i=1;i<=6;i++){reforce[ne][i]=refsave[ne][i];}}
        for(i=1;i<=nt;i++){dist[i]=dsav[i];}
        /*初期処理 {df}=[KT]{dis0}*/
            for(i=1;i<=nt;i++){for(j=1;j<=nt;j++){tsm[i][j]=0.0;}}
            for(ne=1;ne<=NELT;ne++){
                FRAMESM(ne,kakom,x,y,dist,Em,AA,AI,reforce,sm);
                ASMAT(ne,2,3,kakom,sm,tsm);
            }
            for(i=1;i<=nt;i++){reac[i]=0.0;dis0[i]=0.0;work[i]=0.0;}
            for(i=1;i<=mm;i++){
                ia=index[i];work[i]=df[ia];reac[i]=work[i];
                for(j=1;j<=mm;j++){ja=index[j];tsm[i][j]=tsm[ia][ja];}
            }
            switch(jdg){
                case 0: /*Cholesky*/
                    ib=IBAND(mm,tsm,fact0);
                    BAND(mm,ib,tsm);
                    jdg=CHBAND10(mm,ib,tsm);/*三角行列*/
                    CHBAND11(mm,ib,tsm,reac);/*前進消去・後退代入*/
                    break;
                case 1: /*Gauss-Jordan : [tsm]=inverse matrix */
                    jdginv=MATINV(mm,tsm);
                    if(jdginv==1)goto epart;
                    for(i=1;i<=mm;i++){
                        reac[i]=0.0;
                        for(j=1;j<=mm;j++){
                            reac[i]=reac[i]+tsm[i][j]*work[j];
                    }}
                    break;
            }
            for(i=1;i<=mm;i++){ia=index[i];dis0[ia]=reac[i];}
            switch(kprog){
                case 1:dlam=CAL_LAM01(nt,deltaS,Spara,dis0,df);break;
                case 2:dlam=CAL_LAM02(nt,deltaS,coefAB,&Spara,dis0,df);break;
            }
            for(i=1;i<=nt;i++){dis[i]=dlam*dis0[i];}
        /*====================================*/
        /*収束過程*/
        /*====================================*/
        do{
            itaration=itaration+1;
            lam=lam+dlam;
            for(i=1;i<=nt;i++){dist[i]=dist[i]+dis[i];ftvec[i]=0.0;}
            /*内力計算*/
            for(ne=1;ne<=NELT;ne++){
                CALST(ne,kakom,x,y,dis,dist,Em,AA,AI,fevec,reforce);
                ASVEC(ne,2,3,kakom,fevec,ftvec);
            }
            /*不平衡力計算*/
            for(i=1;i<=nt;i++){reac[i]=ff[i]+lam*df[i]-ftvec[i];ftvec[i]=reac[i];}
            for(i=1;i<=mm;i++){ia=index[i];reac[i]=reac[ia];dis[i]=dis[ia];}
            /*収束判定*/
            dtest=0.0;ftest=0.0;
            for(i=1;i<=mm;i++){dtest=dtest+fabs(dis[i]);ftest=ftest+fabs(reac[i]);}
            dtest0=dtest0+dtest;
            ftest0=ftest0+ftest;
            printf("iload=%d itaration=%d dtest=%e ftest=%e lam=%.3e\n",iload,itaration,dtest,ftest,lam);
            /*=====================================================*/
            if(dtest<1e-10)break;
            if(dtest<dtest0*1.0e-4&&ftest<ftest0*1.0e-4)break;
            /*=====================================================*/
            for(i=1;i<=nt;i++){
                dis0[i]=0.0;disr[i]=0.0;reac[i]=0.0;
                for(j=1;j<=nt;j++){tsm[i][j]=0.0;}
            }
            for(ne=1;ne<=NELT;ne++){
                FRAMESM(ne,kakom,x,y,dist,Em,AA,AI,reforce,sm);
                ASMAT(ne,2,3,kakom,sm,tsm);
            }
            for(i=1;i<=mm;i++){
                ia=index[i];reac[i]=df[ia];work[i]=reac[i];
                for(j=1;j<=mm;j++){ja=index[j];tsm[i][j]=tsm[ia][ja];}
            }
            switch(jdg){
                case 0: /*Cholesky*/
                    ib=IBAND(mm,tsm,fact0);
                    BAND(mm,ib,tsm);
                    jdg=CHBAND10(mm,ib,tsm);
                    if(jdg==1)goto again;
                    CHBAND11(mm,ib,tsm,reac);/*前進消去・後退代入*/
                    for(i=1;i<=mm;i++){ia=index[i];dis0[ia]=reac[i];reac[i]=ftvec[ia];}
                    CHBAND11(mm,ib,tsm,reac);/*前進消去・後退代入*/
                    break;
                case 1: /*Gauss-Jordan [tsm]=inverse matrix */
                    jdginv=MATINV(mm,tsm);
                    if(jdginv==1)goto epart;
                    for(i=1;i<=mm;i++){
                        reac[i]=0.0;
                        for(j=1;j<=mm;j++){
                            reac[i]=reac[i]+tsm[i][j]*work[j];
                    }}
                    for(i=1;i<=mm;i++){ia=index[i];dis0[ia]=reac[i];work[i]=ftvec[ia];}
                    for(i=1;i<=mm;i++){
                        reac[i]=0.0;
                        for(j=1;j<=mm;j++){
                            reac[i]=reac[i]+tsm[i][j]*work[j];
                    }}
                    break;
            }
            for(i=1;i<=mm;i++){ia=index[i];disr[ia]=reac[i];}
            dlam=CAL_LAMNEW(nt,Spara,dis0,disr,df);
            for(i=1;i<=nt;i++){dis[i]=dlam*dis0[i]+disr[i];}
        /*====================================*/
        }while(itaration<=maxita);
        /*====================================*/
        if(maxita<itaration)nocon=1;
        /*釣合点確定*/
        for(i=1;i<=nt;i++){ff[i]=ff[i]+lam*df[i];dsav[i]=dist[i];}  /*荷重変位確定*/
        for(ne=1;ne<=NELT;ne++){for(i=1;i<=6;i++){refsave[ne][i]=reforce[ne][i];}}  /*要素内力確定*/
        DATAOUT(fnameW,NODT,NELT,iload,itaration,nocon,deltaS,lam,dtest,ftest,x,y,dist,ff,ftvec,reforce); /*印刷*/
        /*荷重・変位が動かなくなったら計算終了*/
        fend1=0.0;dend1=0.0;
        for(i=1;i<=nt;i++){
            fend1=fend1+fabs(ff[i]);
            dend1=dend1+fabs(dist[i]);
        }
        if(2<=iload){
            if(((1.0-1e-6)*fend0<fend1&&fend1<(1.0+1e-6)*fend0)&&((1.0-1e-6)*fend0<fend1&&fend1<(1.0+1e-6)*fend0)){
                jflagc=1;/*荷重・変位変化停止*/
                goto epart;
            }
        }
        fend0=fend1;dend0=dend1;
        /*弧長変更判定*/
        /*初回の増分荷重に対する増分荷重の比率が1/100未満となったら弧長を coefS 倍にして次回増分荷重を算定*/
        if(iload==1){
            fjdg0=0.0;
            for(i=1;i<=nt;i++){if(fjdg0<fabs(lam*df[i]))fjdg0=fabs(lam*df[i]);}
        }
        fjdg=0.0;
        for(i=1;i<=nt;i++){if(fjdg<fabs(lam*df[i]))fjdg=fabs(lam*df[i]);}
        if(fjdg<fjdg0*0.01)deltaS=coefS*deltaS;
        if(deltaS0*20.0<deltaS)deltaS=deltaS0*20.0;
    /*******************************************************/
    }while(nocon==0&&iload<maxload);
    /*******************************************************/
epart:
    if(nocon==1)jflagc=2;/*未収束*/
    if(jdginv==1)jflagc=3;/*剛性行列特異*/
    /*計算完了*/
    t_end=clock();
    cal_time=(double)(t_end-t_sta)/CLOCKS_PER_SEC;
    time(&jikoku);
    lt=localtime(&jikoku);
    p=asctime(lt);

    fout=fopen(fnameW,"a+");
    switch(jflagc){
        case 0:fprintf(fout,"jflagc=%d : Normal completion\n",jflagc);break;
        case 1:fprintf(fout,"jflagc=%d : Stand still\n"      ,jflagc);break;
        case 2:fprintf(fout,"jflagc=%d : No convergence\n"   ,jflagc);break;
        case 3:fprintf(fout,"jflagc=%d : Singular matrix\n"  ,jflagc);break;
    }
    fprintf(fout,"NODT=%d nt=%d mm=%d ib=%d\n",NODT,nt,mm,ib);
    fprintf(fout,"deltaS=%g jdg=%d iload=%d itaration=%d\n",deltaS,jdg,iload,itaration);
    fprintf(fout,"Calculation time=%.3f(sec)\n",cal_time);
    fprintf(fout,"%s\n",p);
    fclose(fout);
 
    switch(jflagc){
        case 0:printf("jflagc=%d : Normal completion\n",jflagc);break;
        case 1:printf("jflagc=%d : Stand still\n"      ,jflagc);break;
        case 2:printf("jflagc=%d : No convergence\n"   ,jflagc);break;
        case 3:printf("jflagc=%d : Singular matrix\n"  ,jflagc);break;
    }
    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;
 
}
/*--------------------------------------------------------------------------------*/
double CAL_LAM01(int nt,double deltaS,double Spara,double dis0[NTmax+1],double df[NTmax+1])
{
    int i;
    double sum1,sum2;

    sum1=0.0;sum2=0.0;
    for(i=1;i<=nt;i++){
        sum1=sum1+dis0[i]*dis0[i];
        sum2=sum2+df[i]*df[i];
    }
    return sqrt(deltaS*deltaS/(sum1+(Spara)*(Spara)*sum2));

}
/*--------------------------------------------------------------------------------*/
double CAL_LAM02(int nt,double deltaS,double coefAB,double *Spara,double dis0[NTmax+1],double df[NTmax+1])
{
    int i;
    double sum1,sum2,sum3;

    sum1=0.0;sum2=0.0;sum3=0.0;
    for(i=1;i<=nt;i++){
        sum1=sum1+dis0[i]*dis0[i];
        sum2=sum2+df[i]*df[i];
        sum3=sum3+dis0[i]*df[i];
    }
    *Spara=sum3/sum2/coefAB;
    return coefAB*sqrt(deltaS*deltaS/(sum1+(*Spara)*(*Spara)*sum2));
}
/*--------------------------------------------------------------------------------*/
double CAL_LAMNEW(int nt,double Spara,double dis0[NTmax+1],double disr[NTmax+1],double df[NTmax+1])
{
    int i;
    double sum1,sum2,sum3;

    sum1=0.0;
    sum2=0.0;
    sum3=0.0;
    for(i=1;i<=nt;i++){
        sum1=sum1+dis0[i]*disr[i];
        sum2=sum2+dis0[i]*dis0[i];
        sum3=sum3+df[i]*df[i];
    }
    return -sum1/(sum2+Spara*Spara*sum3);
}
/*--------------------------------------------------------------------------------*/
void DATAOUT(char fnameW[50],int NODT,int NELT,
    int iload,int itaration,int nocon,double deltaS,double lam,double dtest,double ftest,
    double x[NDmax+1],double y[NDmax+1],double dist[NTmax+1],
    double ff[NTmax+1],double ftvec[NTmax+1],double reforce[NEmax+1][7])
{
    FILE *fout;
    int i,ne;
    char dat[256]="";
    char sv[50]="";

    fout=fopen(fnameW,"a+");
    switch(nocon){
        case 0:fprintf(fout,"***** Output *****\n");break;
        case 1:fprintf(fout,"##### Output #####\n");break;
    }
    fprintf(fout,"iload,itaration,nocon,deltaS,lam,dtest,ftest\n");
        sprintf(sv,"%d",iload)    ;strcpy(dat,sv);
        sprintf(sv,"%d",itaration);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%d",nocon)    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",deltaS)      ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",lam)      ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",dtest)    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ftest)    ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    fprintf(fout,"*Node values\n");
    fprintf(fout,"node,x-coord,y-coord,x-dir,y-dir,z-dir,ff-x,ff-y,ff-z,ftvec-x,ftvec-y,ftvec-z\n");
    for(i=1;i<=NODT;i++){
        sprintf(sv,"%d",i)               ;strcpy(dat,sv);
        sprintf(sv,"%e",x[i]+dist[3*i-2]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",y[i]+dist[3*i-1]);strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",dist[3*i-2])     ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",dist[3*i-1])     ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",dist[3*i])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ff[3*i-2])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ff[3*i-1])       ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ff[3*i])         ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ftvec[3*i-2])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ftvec[3*i-1])    ;strcat(dat,",");strcat(dat,sv);
        sprintf(sv,"%e",ftvec[3*i])      ;strcat(dat,",");strcat(dat,sv);
        fprintf(fout,"%s\n",dat);
    }
    fprintf(fout,"*Stress-resultant\n");
    fprintf(fout,"element,Ni,Si,Mi,Nj,Sj,Mj\n");
    for(ne=1;ne<=NELT;ne++){
        sprintf(sv,"%d",ne);strcpy(dat,sv);
        for(i=1;i<=6;i++){
            sprintf(sv,"%e",reforce[ne][i]);strcat(dat,",");strcat(dat,sv);
        }
        fprintf(fout,"%s\n",dat);
    }
    fclose(fout);
}
/*--------------------------------------------------------------------------------*/
void FRAMESM(int ne,int kakom[NEmax+1][3],double x[NDmax+1],double y[NDmax+1],double dist[NTmax+1],
    double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double reforce[NEmax+1][7],double sm[7][7])
{
    int i,j,k;
    double hen[7][7];
    double sme[7][7];
    double smv[7][7];
    double w[7][7];
    double AL,P;
    double xx1,yy1,xx2,yy2;

    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            smv[i][j]=0.0;
            sm[i][j]=0.0;
        }
    }
    /*軸力計算*/
    i=kakom[ne][1];
    j=kakom[ne][2];
    xx1=x[i]+dist[3*i-2];
    yy1=y[i]+dist[3*i-1];
    xx2=x[j]+dist[3*j-2];
    yy2=y[j]+dist[3*j-1];
    AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1));
    P=0.5*(reforce[ne][4]-reforce[ne][1]);

    /*幾何剛性行列*/
    smv[1][1]=0.0;
    smv[2][1]=0.0;
    smv[2][2]=6.0/5/AL;
    smv[3][1]=0.0;
    smv[3][2]=1.0/10.0;
    smv[3][3]=2.0*AL/15.0;
    smv[4][1]=0.0;
    smv[4][2]=0.0;
    smv[4][3]=0.0;
    smv[4][4]=0.0;
    smv[5][1]=0.0;
    smv[5][2]=-6.0/5.0/AL;
    smv[5][3]=-1.0/10.0;
    smv[5][4]=0.0;
    smv[5][5]=6.0/5.0/AL;
    smv[6][1]=0.0;
    smv[6][2]=1.0/10.0;
    smv[6][3]=-AL/30.0;
    smv[6][4]=0.0;
    smv[6][5]=-1.0/10.0;
    smv[6][6]=2.0*AL/15.0;
    for(i=1;i<=6;i++){
        for(j=i;j<=6;j++){
            smv[i][j]=smv[j][i];
        }
    }
    /*線形弾性剛性行列*/
    ELMSM(ne,xx1,yy1,xx2,yy2,Em,AA,AI,sme);
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            sm[i][j]=P*smv[i][j]+sme[i][j];
        }
    }
    HENMX(ne,xx1,yy1,xx2,yy2,hen);
    /*[w]=[T]/**[k]*/
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            w[i][j]=0.0;
            for(k=1;k<=6;k++){
                w[i][j]=w[i][j]+hen[k][i]*sm[k][j];
            }
        }
    }
    /*[w]*[T]=[T]/*[k][T]：[sm]は全体座標系での要素剛性行列*/
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            sm[i][j]=0.0;
            for(k=1;k<=6;k++){
                sm[i][j]=sm[i][j]+w[i][k]*hen[k][j];
            }
        }
    }

}
/*--------------------------------------------------------------------------------*/
void ELMSM(int ne,double xx1,double yy1,double xx2,double yy2,
    double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double sme[7][7])
{
    /*要素剛性行列*/
    int i,j;
    double AL,EA,EI;

    AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1));
    EA=Em[ne]*AA[ne];
    EI=Em[ne]*AI[ne];
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            sme[i][j]=0.0;
        }
    }
    sme[1][1]=EA/AL;
    sme[2][2]=12.0*EI/(AL*AL*AL);
    sme[3][2]=6.0*EI/(AL*AL);
    sme[3][3]=4.0*EI/AL;
    sme[4][1]=-EA/AL;
    sme[4][4]=sme[1][1];
    sme[5][2]=-sme[2][2];
    sme[5][3]=-sme[3][2];
    sme[5][5]=sme[2][2];
    sme[6][2]=sme[3][2];
    sme[6][3]=2.0*EI/AL;
    sme[6][5]=-sme[3][2];
    sme[6][6]=sme[3][3];
    /*対称項セット*/
    for(i=1;i<=6;i++){
        for(j=1;j<=i;j++){
            sme[j][i]=sme[i][j];
        }
    }
}
/*--------------------------------------------------------------------------------*/
void CALST(int ne, int kakom[NEmax+1][3], double x[NDmax+1],double y[NDmax+1],double dis[NTmax+1],
    double dist[NTmax+1],double Em[NEmax+1],double AA[NEmax+1],double AI[NEmax+1],double fevec[7],double reforce[NEmax+1][7])
{
    /*断面力算定*/
    int i,k,ii,jj,ia;
    double hen[7][7],sme[7][7];
    double AL,BL,deltai,deltaj,wd[7],wde[7],ws[7];
    double xx1,yy1,xx2,yy2;

    ROTATE(ne,kakom,x,y,dis,dist,&deltai,&deltaj);

    /*要素節点変位増分（全体座標系）*/
    for(i=1;i<=2;i++){
        ia=kakom[ne][i];
        wd[3*i-2]= dis[3*ia-2];
        wd[3*i-1]= dis[3*ia-1];
        wd[3*i]  = dis[3*ia];
    }
    ii=kakom[ne][1];
    jj=kakom[ne][2];
    /*前回要素座標系からの変位増分*/
    xx1=x[ii]+dist[3*ii-2]-dis[3*ii-2];
    yy1=y[ii]+dist[3*ii-1]-dis[3*ii-1];
    xx2=x[jj]+dist[3*jj-2]-dis[3*jj-2];
    yy2=y[jj]+dist[3*jj-1]-dis[3*jj-1];
    HENMX(ne,xx1,yy1,xx2,yy2,hen);
    for(i=1;i<=6;i++){
        wde[i]=0.0;
        for(k=1;k<=6;k++){
            wde[i]=wde[i]+hen[i][k]*wd[k];
    }}
    /*前回部材長*/
    xx1=x[ii]+dist[3*ii-2]-dis[3*ii-2];
    yy1=y[ii]+dist[3*ii-1]-dis[3*ii-1];
    xx2=x[jj]+dist[3*jj-2]-dis[3*jj-2];
    yy2=y[jj]+dist[3*jj-1]-dis[3*jj-1];
    BL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1));
    /*今回部材長*/
    xx1=x[ii]+dist[3*ii-2];
    yy1=y[ii]+dist[3*ii-1];
    xx2=x[jj]+dist[3*jj-2];
    yy2=y[jj]+dist[3*jj-1];
    AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1));
    wd[1]=0.0;
    wd[2]=0.0;
    wd[3]=deltai;
    wd[4]=AL-BL;
    wd[5]=0.0;
    wd[6]=deltaj;
    /*要素内力増分(前回微小変形要素剛性使用)*/
    xx1=x[ii]+dist[3*ii-2]-dis[3*ii-2];
    yy1=y[ii]+dist[3*ii-1]-dis[3*ii-1];
    xx2=x[jj]+dist[3*jj-2]-dis[3*jj-2];
    yy2=y[jj]+dist[3*jj-1]-dis[3*jj-1];
    ELMSM(ne,xx1,yy1,xx2,yy2,Em,AA,AI,sme);
    for(i=1;i<=6;i++){
        ws[i]=0.0;
        for(k=1;k<=6;k++){
            ws[i]=ws[i]+sme[i][k]*wd[k];
    }}
    /*断面力保存*/
    for(i=1;i<=6;i++){
        reforce[ne][i]=reforce[ne][i]+ws[i];
    }
    /*全体座標系での変形後要素内力*/
    xx1=x[ii]+dist[3*ii-2];
    yy1=y[ii]+dist[3*ii-1];
    xx2=x[jj]+dist[3*jj-2];
    yy2=y[jj]+dist[3*jj-1];
    HENMX(ne,xx1,yy1,xx2,yy2,hen);
    for(i=1;i<= 6;i++){
        fevec[i]=0.0;
        for(k=1;k<=6;k++){
            fevec[i]=fevec[i]+hen[k][i]*reforce[ne][k];
    }}
}
/*--------------------------------------------------------------------------------*/
void ROTATE(int ne, int kakom[NEmax+1][3], double x[NDmax+1],double y[NDmax+1],
    double dis[NTmax+1],double dist[NTmax+1],double *deltai,double *deltaj)
{
    /*剛体回転除去*/
    int i,j,k,ii,jj,ia;
    double hen[7][7];
    double wd[7],wde[7];
    double xx1,yy1,xx2,yy2;
    double BL,Ti,Tj,theta1i,theta1j,theta2i,theta2j;

    ii=kakom[ne][1];
    jj=kakom[ne][2];
    /*前回における初期部材座標系からの変位量*/
    for(i=1;i<=2;i++){
        ia=kakom[ne][i];
        wd[3*i-2]=dist[3*ia-2]-dis[3*ia-2];
        wd[3*i-1]=dist[3*ia-1]-dis[3*ia-1];
        wd[3*i]  =dist[3*ia]  -dis[3*ia];
    }
    /*初期部材長さ*/
    xx1=x[ii];
    yy1=y[ii];
    xx2=x[jj];
    yy2=y[jj];
    BL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1));
    HENMX(ne,xx1,yy1,xx2,yy2,hen);
    for(i=1;i<=6;i++){
        wde[i]=0.0;
        for(k=1;k<=6;k++){
            wde[i]=wde[i]+hen[i][k]*wd[k];
    }}
    /*前回の内力計算用回転角*/
    Ti=tan(wde[3]);
    Tj=tan(wde[6]);
    theta1i=((BL+wde[4]-wde[1])*Ti-(wde[5]-wde[2]))/((BL+wde[4]-wde[1])+(wde[5]-wde[2])*Ti);
    theta1j=((BL+wde[4]-wde[1])*Tj-(wde[5]-wde[2]))/((BL+wde[4]-wde[1])+(wde[5]-wde[2])*Tj);

    /*今回における初期部材座標系からの変位量*/
    for(i=1;i<=2;i++){
        ia=kakom[ne][i];
        wd[3*i-2]=dist[3*ia-2];
        wd[3*i-1]=dist[3*ia-1];
        wd[3*i]  =dist[3*ia]  ;
    }
    for(i=1;i<=6;i++){
        wde[i]=0.0;
        for(k=1;k<=6;k++){
            wde[i]=wde[i]+hen[i][k]*wd[k];
    }}
    /*前回の内力計算用回転角*/
    Ti=tan(wde[3]);
    Tj=tan(wde[6]);
    theta2i=((BL+wde[4]-wde[1])*Ti-(wde[5]-wde[2]))/((BL+wde[4]-wde[1])+(wde[5]-wde[2])*Ti);
    theta2j=((BL+wde[4]-wde[1])*Tj-(wde[5]-wde[2]))/((BL+wde[4]-wde[1])+(wde[5]-wde[2])*Tj);

    /*内力計算用回転角増分*/
    *deltai=theta2i-theta1i;
    *deltaj=theta2j-theta1j;

}
/*--------------------------------------------------------------------------------*/
void HENMX(int ne,double xx1,double yy1,double xx2,double yy2,double hen[7][7])
{
    /*座標変換行列*/
    int i,j;
    double AL,cs,sn;

    AL=sqrt((xx2-xx1)*(xx2-xx1)+(yy2-yy1)*(yy2-yy1));
    cs=(xx2-xx1)/AL;
    sn=(yy2-yy1)/AL;
    for(i=1;i<=6;i++){
        for(j=1;j<=6;j++){
            hen[i][j]=0.0;
        }
    }
    hen[1][1]=cs;
    hen[1][2]=sn;
    hen[2][1]=-sn;
    hen[2][2]=cs;
    hen[3][3]=1.0;
    for(i=1;i<=3;i++){
        for(j=1;j<=3;j++){
            hen[i+3][j+3]=hen[i][j];
        }
    }
}
/*--------------------------------------------------------------------------------*/
void ASMAT(int ne,int nod,int nhen,int kakom[NEmax+1][3],double sm[7][7],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 ASVEC(int ne,int nod,int nhen,int kakom[NEmax+1][3],double fevec[7],double ftvec[NTmax+1])
{
    /*荷重ベクトル組立[nod：1要素節点数，nhen：1節点の自由度]*/
    int i,k,ki,ks,kik,ksk;

    for(i=1;i<=nod;i++){
        ki=(kakom[ne][i]-1)*nhen;
        ks=(i-1)*nhen;
        for(k=1;k<=nhen;k++){
            kik=ki+k;
            ksk=ks+k;
            ftvec[kik]=ftvec[kik]+fevec[ksk];
        }
    }

}
/*--------------------------------------------------------------------------------*/
int 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{
                    if(z<0)return 1;
                    array[i][1]=sqrt(z);
                }
            }else{
                array[i][j]=z/array[i][1];
            }
        }
    }
    return 0;
}
/*--------------------------------------------------------------------------------*/
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];
        }
    }
}
/*--------------------------------------------------------------------------------*/
int MATINV(int n,double ary[NTmax+1][NTmax+1])
{
    /*****************************************/
    /*  逆行列計算                           */
    /*  山田嘉昭著：マトリックス法材料力学   */
    //****************************************/
    int i,j,k,l,l1,irow=0,icol=0,jrow=0,jcol=0;
    double amax,s,t;
    double determ=1.0;
    double factor=1.0e-30;
    int ipivot[n+1];
    double pivot[n+1];
    int index[n+1][3];

    for(j=1;j<=n;j++){
        ipivot[j]=0;
    }

    for(i=1;i<=n;i++){

        amax=0.0;
        for(j=1;j<=n;j++){
            if(ipivot[j]-1!=0){
                for(k=1;k<=n;k++){
                    if(0<ipivot[k]-1)return 0;
                    if(ipivot[k]-1<0){
                        if(fabs(amax)-fabs(ary[j][k])<0){
                            irow=j;
                            icol=k;
                            amax=ary[j][k];
                        }
                    }
                }
            }
        }
        ipivot[icol]=ipivot[icol]+1;

        if((irow+icol)%2!=0){
            determ=-determ;
        }
        if(irow-icol!=0){
            for(l=1;l<=n;l++){
                s=ary[irow][l];
                ary[irow][l]=ary[icol][l];
                ary[icol][l]=s;
            }
        }
        index[i][1]=irow;
        index[i][2]=icol;
        pivot[i]=ary[icol][icol];
        if(fabs(pivot[i])<factor)return 1;
        if(determ<=1.0E+30){
            determ=determ*pivot[i];
        }
        ary[icol][icol]=1.0;
        for(l=1;l<=n;l++){
            ary[icol][l]=ary[icol][l]/pivot[i];
        }

        for(l1=1;l1<=n;l1++){
            if(l1-icol!=0){
                t=ary[l1][icol];
                ary[l1][icol]=0.0;
                for(l=1;l<=n;l++){
                    ary[l1][l]=ary[l1][l]-ary[icol][l]*t;
                }
            }
        }

    }

    for(i=1;i<=n;i++){
        l=n+1-i;
        if(index[l][1]-index[l][2]!=0){
            jrow=index[l][1];
            jcol=index[l][2];
            for(k=1;k<=n;k++){
                s=ary[k][jrow];
                ary[k][jrow]=ary[k][jcol];
                ary[k][jcol]=s;
            }
        }
    }

    return 0;

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