Examples of QRdcmp

com.nr.la.QRdcmp
QR Decomposition - PTC There is another matrix factorization that is sometimes very useful, the so-called QR decomposition, A=Q*R Here R is upper triangular, while Q is orthogonal, that is, QT*Q = 1 QR decomposition can be used to solve systems of linear equations. To solve A*x=b first form QT*b and then solve R*x=QT*b by backsubstitution. Copyright (C) Numerical Recipes Software 1986-2007 Java translation Copyright (C) Huang Wen Hui 2012 @author hwh

Examples of com.nr.la.QRdcmp


    


    // Test QRdcmp
    System.out.println("Testing QRdcmp");
    QRdcmp aqr = new QRdcmp(a);
    sbeps = 5.e-14;
    localflag = maxel(matsub(matmul(transpose(aqr.qt),aqr.r),a)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** QRdcmp: QR decomposition was unsuccessful");
      
    }


    localflag = maxel(matsub(matmul(transpose(aqr.qt),aqr.qt),ident(a.length,1.))) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** QRdcmp: Matrix aqr.qt is not orthogonal");
      
    }


    ranvec(r);
    aqr.solve(r,y);
    localflag = maxel(vecsub(matmul(a,y),r)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** QRdcmp: Error in solve() method");
      
    }


    // test update
    double[][] aa = buildMatrix(a);
    double[] s = new double[aa.length],t = new double[aa.length];
    ranvec(s); ranvec(t);
    double[] u = new double[aa.length],v = buildVector(t);
    aqr.qtmult(s,u);
    aqr.update(u,v);
    for (i=0;i<aa.length;i++) for (j=0;j<aa.length;j++) a[i][j] = aa[i][j]+s[i]*t[j];
    localflag = maxel(matsub(matmul(transpose(aqr.qt),aqr.r),a)) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** QRdcmp: Error in method qtmult() or update()");

View Full Code Here

Examples of com.nr.la.QRdcmp

        1.0,-1.0,  5.0};
    double[][] aa = NRUtil.buildMatrix(3,3, a);
    double[] b = new double[]{2.0, 4.0,1.0};
    double[] x = new double[3];


    QRdcmp qr = new QRdcmp(aa);
    qr.solve(b, x);
    System.out.println(NRUtil.toString(x));
  }

View Full Code Here

Examples of com.nr.la.QRdcmp

    boolean restrt,skip;
    int i,its,j,n=x.length;
    double den,fold,stpmax,sum,temp,test;
    doubleW f = new doubleW(0);
    double[] fvcold =new double[n],g=new double[n],p=new double[n],s=new double[n],t=new double[n],w=new double[n],xold=new double[n];
    QRdcmp qr=null;
    NRfmin fmin = new NRfmin(vecfunc);
    NRfdjac fdjac = new NRfdjac(vecfunc);
    f.val=fmin.funk(x);  // XXX init fmin.fvec first.
    double[] fvec=fmin.fvec;
    test=0.0;
    for (i=0;i<n;i++)
      if (abs(fvec[i]) > test) test=abs(fvec[i]);
    if (test < 0.01*TOLF) {
      check.val=false;
      return;
    }
    for (sum=0.0,i=0;i<n;i++) sum += SQR(x[i]);
    stpmax=STPMX*max(sqrt(sum),n);
    restrt=true;
    for (its=1;its<=MAXITS;its++) {
      if (restrt) {
        qr=new QRdcmp(fdjac.get(x,fvec));
        if (qr.sing) {
          double[][] one = new double[n][n];
          for (i=0;i<n;i++) one[i][i]=1.0;
          qr=new QRdcmp(one);
        }
      } else {
        for (i=0;i<n;i++) s[i]=x[i]-xold[i];
        for (i=0;i<n;i++) {
          for (sum=0.0,j=i;j<n;j++) sum += qr.r[i][j]*s[j];
          t[i]=sum;
        }
        skip=true;
        for (i=0;i<n;i++) {
          for (sum=0.0,j=0;j<n;j++) sum += qr.qt[j][i]*t[j];
          w[i]=fvec[i]-fvcold[i]-sum;
          if (abs(w[i]) >= EPS*(abs(fvec[i])+abs(fvcold[i]))) skip=false;
          else w[i]=0.0;
        }
        if (!skip) {
          qr.qtmult(w,t);
          for (den=0.0,i=0;i<n;i++) den += SQR(s[i]);
          for (i=0;i<n;i++) s[i] /= den;
          qr.update(t,s);
          if (qr.sing) throw new IllegalArgumentException("singular update in broydn");
        }
      }
      qr.qtmult(fvec,p);
      for (i=0;i<n;i++)
        p[i] = -p[i];
      for (i=n-1;i>=0;i--) {
        for (sum=0.0,j=0;j<=i;j++) sum -= qr.r[j][i]*p[j];
        g[i]=sum;
      }
      for (i=0;i<n;i++) {
        xold[i]=x[i];
        fvcold[i]=fvec[i];
      }
      fold=f.val;
      qr.rsolve(p,p);
      double slope=0.0;
      for (i=0;i<n;i++) slope += g[i]*p[i];
      if (slope >= 0.0) {
        restrt=true;
        continue;

View Full Code Here

TOP

All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.