normdd_s.c

void edit_norm_dd(icon,idN0,id_slips,N0t,slips_t,am,lo,idlo,idN00,N00,N0tb,
                  m00,ici,c_dd,an,b1,an_s,b1_s,ltpl,ltpl_s,
                  an_s1,b1_s1,ltpl_s1,an_k,b1_k,ltpl_k)
int icon[I10*Isat],ici[Iunknown];
long idN0[Iunknown],id_slips[Islips],N0t[Iunknown];
double slips_t[Islips],am[Isato*Ista*Iunknown],lo[Isato*Ista];
long idlo[Iunknown],idN00[Iunknown],N0tb[Iunknown];
double N00[Iunknown],m00;
double c_dd[Nlow],an[Nlow],b1[N],an_s[Nlow],b1_s[N],*ltpl,*ltpl_s;
double an_k[Nlow],b1_k[N],an_s1[Nlow],b1_s1[N],*ltpl_k,*ltpl_s1;
/* 函數(shù) edit_norm 應(yīng)用從 eq_dd_s.c 輸入的矩陣重新排列觀測方程 */ 
{
     long i,j,j1,j2,k,k1;
     double sum,sum1;
/* 刪掉未知數(shù)和應(yīng)用已知數(shù)解并組成 idN00，用3 sigm刪除obs.，為dd 數(shù)據(jù)設(shè)權(quán)矩陣 */
     j1 = 2*(icon[200] + 1);
     for (i=0;i<icon[200];i++) 
      { for (j=0;j<=i;j++) 
         { c_dd[j+i*(i+1)/2] = -1./j1;
           if(j==i) c_dd[j+i*(i+1)/2] = ((double)icon[200])/j1;}}
/* 初始化 an[], b1[], *ltpl, an_k[], b1_k[], *ltpl_k */
     for (i=0;i<Iunknown;i++) 
      {for(j=0;j<=i;j++) {an[j+i*(i+1)/2]=0.;an_k[j+i*(i+1)/2] = 0.;}
        b1[i] = 0.; b1_k[i]=0.;}
       *ltpl = 0.; *ltpl_k = 0.;
/* 組成法方程， dd 數(shù)據(jù)形式，計算 N = amT*P*am和 B1 = amT*P*lo */
     j1 = 2*(icon[200] + 1);
     for (i=0;i<icon[201];i++) 
      { for (j=0;j<=i;j++) 
         { sum = 0.; 
           for (k=0;k<icon[200];k++) 
            { for (k1=0;k1<icon[200];k1++) 
              { j2 = -1; if(k==k1) j2 = icon[200];
                sum += am[i+k*Iunknown]*j2*am[j+k1*Iunknown]; }} 
           an[j+i*(i+1)/2] = sum/j1; }
        sum1 = 0.; 
        for (j=0;j<icon[200];j++) 
     { for (k=0;k<icon[200];k++) 
       {   j2 = -1; if(j==k) j2 = icon[200];
           sum1 += am[i+j*Iunknown]*j2*lo[k]; }}
        b1[i] = sum1/j1; }              
     sum = 0.;
     for (j=0;j<icon[200];j++) 
     { for (k=0;k<icon[200];k++) 
       { j2 = -1; if(j==k) j2 = icon[200];
         sum += lo[j]*j2*lo[k]; }}
     *ltpl = sum/j1;
        for (i=0;i<icon[201];i++) 
        { for (j=0;j<=i;j++) 
          {   an_s[j+i*(i+1)/2]+=an[j+i*(i+1)/2];
              an_s1[j+i*(i+1)/2]=an_s[j+i*(i+1)/2]; }
           b1_s[i]+=b1[i]; b1_s1[i]=b1_s[i]; }              
        *ltpl_s += *ltpl; *ltpl_s1 = *ltpl_s; icon[202]+=icon[200];
        if(icon[202]>=10000){icon[202]-=10000; icon[207]+=1;}
/* 為控制輸出數(shù)據(jù) */
}