惊叹

freddie_1999

这些天在整个类似横观各项同性本构的2D开发，惊讶于ITASCA程序开发之高超，聊聊数语即能解决繁琐的矩阵求逆过程，我在BAIDU上搜了好久，提供的源程序要不就是错的，要不就是繁琐之极。不免敬意油然而生，故立此贴以示崇拜！！！！
bool Mymodel::xmatinv(double b[6][6]) {
  bool flag;
  double bii,bji;
  int n=6;
  int i,j,k;
  flag = true;
  for (i=0; i<n; i++) {
    bii = b[i][i];
    if (bii == 0.0) {
      flag = false;
      return(flag);
    }
    for (k=0; k<n; k++) b[i][k] = b[i][k] / bii;
    b[i][i] = 1.0 / bii;
    for (j=0; j<n; j++) {
      if (j != i) {
        bji = b[j][i];
        for (k=0; k<n; k++) b[j][k] = b[j][k] - bji * b[i][k];
        b[j][i] = - bji * b[i][i];
      }
    }
  }
  return(flag);
}

[[i] 本帖最后由 benjackxu 于 2006-12-21 11:28 编辑 [/i]]

sch

呵呵！这是C++的功效！看看我的一段雅克比求矩阵特征值代码：
//雅可比求特征向量和特征值
int eejcb(double **a,//求解矩阵
                  int n,//矩阵行列数n*n
                  double **v,//返回的特征向量矩阵
                  double eps,//收敛精度
                  int jt//迭带次数
                  )
  {
          int i,j,p,q,l;
          double fm,cn,sn,omega,x,y,d;
          l=1;
                for (i=0; i<=n-1; i++)
                {
                        v=1.0;
                        for (j=0; j<=n-1; j++)
                                if (i!=j)
                                        v[j]=0.0;
                }
                while (1==1)
                {
                        fm=0.0;
                        for (i=0; i<=n-1; i++)
                                for (j=0; j<=n-1; j++)
                                {
                                        d=fabs(a[j]);
                                        if ((i!=j)&&(d>fm))
                                                {
                                                fm=d;
                                                p=i;
                                                q=j;
                                                }
                                }
                        if (fm<eps)
                                return(1);
                        if (l>jt)
                                return(-1);
                        l=l+1;
                        //u=p*n+q; w=p*n+p; t=q*n+p; s=q*n+q;
                       
                        x=-a[p][q]; y=(a[q][q]-a[p][p])/2.0;
                        omega=x/sqrt(x*x+y*y);
                        if (y<0.0)
                                omega=-omega;
                        sn=1.0+sqrt(1.0-omega*omega);
                        sn=omega/sqrt(2.0*sn);
                        cn=sqrt(1.0-sn*sn);
                        fm=a[p][p];
                        a[p][p]=fm*cn*cn+a[q][q]*sn*sn+a[p][q]*omega;
                        a[q][q]=fm*sn*sn+a[q][q]*cn*cn-a[p][q]*omega;
                        a[p][q]=0.0;
                        a[q][p]=0.0;
                       
                 for (j=0; j<=n-1; j++)
                        if ((j!=p)&&(j!=q))
                        {
                //                u=p*n+j; w=q*n+j;
                                fm=a[p][j];
                                a[p][j]=fm*cn+a[q][j]*sn;
                                a[q][j]=-fm*sn+a[q][j]*cn;
                         }
                for (i=0; i<=n-1; i++)
                        if ((i!=p)&&(i!=q))
            {
                        //        u=i*n+p; w=i*n+q;
                                fm=a[p];
                                a[p]=fm*cn+a[q]*sn;
                                a[q]=-fm*sn+a[q]*cn;
            }
        for (i=0; i<=n-1; i++)
          {
                //        u=i*n+p; w=i*n+q;
            fm=v[p];
            v[p]=fm*cn+v[q]*sn;
            v[q]=-fm*sn+v[q]*cn;
          }
               
      }//while
               
    return(1);
       
}

lakewater

呵呵

benjackxu

加分加分加分！！！！崇拜死了，俺也学学C++。

[ 本帖最后由 benjackxu 于 2006-12-21 11:27 编辑 ]