Changeset 2050 for branches/newmole/source

Timestamp:

05/12/08 16:59:54 (8 months ago)

Author:

rjrw

Message:

Further tests on:

the consistency of the Householder matrix with the Jacobian in
mole_solve

underlying accuracy of projection in and out of the Householder basis
in TestHouseholder?

Location:

branches/newmole/source

Files:

• mole_solve.cpp (modified) (3 diffs)
• TestHouseholder.cpp (modified) (2 diffs)

branches/newmole/source/mole_solve.cpp

@@ -485 +485 @@
                                 }
                         }
+
                         valarray<double> b2(mole.num_calc);
                         for(i=0;i<mole.num_calc;i++)
@@ -493 +494 @@
                         {
                                 valarray<double> v(mole.num_calc-nelem),
-                                        c(mole.num_calc-nelem),
+                                        b3(mole.num_calc-nelem),
                                         p(mole.num_calc-nelem);
                                 for(i=0;i<mole.num_calc-nelem;i++)
                                 {
                                         v[i] = u[i+nelem][nelem];
-                                        c[i] = b2[i+nelem];
-                                }
-                                p = project(c,v,betah[nelem]);
+                                        b3[i] = b2[i+nelem];
+                                }
+                                p = project(b3,v,betah[nelem]);
                                 for (i=0;i<mole.num_calc-nelem;i++)
                                 {
@@ -513 +514 @@
                         {
                                 fprintf(ioQQQ,"Remain %ld proj %g\n",nelem,b2[nelem]);
+                        }
+
+                        for(long j=0;j<mole.num_calc;++j)
+                        {
+                                for(i=0;i<mole.num_calc;++i)
+                                {
+                                        b2[i] = c[j][i];
+                                }
+                                for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                                {
+                                        valarray<double> v(mole.num_calc-nelem),
+                                                b3(mole.num_calc-nelem),
+                                                p(mole.num_calc-nelem);
+                                        for(i=0;i<mole.num_calc-nelem;i++)
+                                        {
+                                                v[i] = u[i+nelem][nelem];
+                                                b3[i] = b2[i+nelem];
+                                        }
+                                        p = project(b3,v,betah[nelem]);
+                                        for (i=0;i<mole.num_calc-nelem;i++)
+                                        {
+                                                b2[i+nelem] = p[i];
+                                        }
+                                }
+                                fprintf(ioQQQ,"Eqn... %ld\n",j);
+                                for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                                {
+                                        fprintf(ioQQQ,"Eqqq %ld proj %g\n",nelem,b2[nelem]);
+                                }
+                                for (;nelem<mole.num_calc;nelem++)
+                                {
+                                        fprintf(ioQQQ,"Remn %ld proj %g\n",nelem,b2[nelem]);
+                                }
                         }
                 }

branches/newmole/source/TestHouseholder.cpp

@@ -43 +43 @@
         TEST_FIXTURE(HouseholderFixture,TestHouseholderOther)
         {
-                valarray<double> v(x.size()),r(x.size());
+                valarray<double> v(x.size()),r(x.size()),b(x.size());
                 double beta;
 
@@ -50 +50 @@
                 CHECK_CLOSE((y*y).sum(),(r*r).sum(),1e-6);
                 CHECK_CLOSE((y*x).sum()/sqrt((x*x).sum()),r[0],1e-6);
+                b = project(r,v,beta);
+                CHECK_CLOSE((b*y).sum(),(y*y).sum(),1e-6);
         }
         TEST_FIXTURE(HouseholderFixture,TestHouseholderProj2)