@@ -754,6 +754,15 @@ int RunReconciliation(DATA* data, threadData_t *threadData, inputData x, matrixD
754754 solveMatrixMultiplication (tmpmatrixC,jacFt.data ,jacF.rows ,Sx.column ,jacFt.rows ,jacFt.column ,tmpmatrixD);
755755 // printMatrix(tmpmatrixD,jacF.rows,jacFt.column,"F*Sx*Ft");
756756 // printMatrix(setc,nsetcvars,1,"c(x,y)");
757+ /*
758+ * Copy tmpmatrixC and tmpmatrixD to avoid loss of data
759+ * when calculating F*
760+ */
761+ matrixData cpytmpmatrixC = {jacF.rows ,Sx.column ,tmpmatrixC};
762+ matrixData cpytmpmatrixD = {jacF.rows ,jacFt.column ,tmpmatrixD};
763+ matrixData tmpmatrixC1 = copyMatrix (cpytmpmatrixC);
764+ matrixData tmpmatrixD1 = copyMatrix (cpytmpmatrixD);
765+
757766 if (ACTIVE_STREAM (LOG_JAC))
758767 {
759768 cout << " Calculations of Matrix (F*Sx*Ft) f* = c(x,y) " << " \n " ;
@@ -784,25 +793,27 @@ int RunReconciliation(DATA* data, threadData_t *threadData, inputData x, matrixD
784793 {
785794 cout << " Calculations of Matrix (F*Sx*Ft) F* = F*Sx " << " \n " ;
786795 cout << " ===============================================" ;
787- printMatrix (tmpmatrixC,jacF.rows ,Sx.column ," F*Sx" );
788- printMatrix (tmpmatrixD,jacF.rows ,jacFt.column ," F*Sx*Ft" );
796+ // printMatrix(tmpmatrixC,jacF.rows,Sx.column,"F*Sx");
797+ // printMatrix(tmpmatrixD,jacF.rows,jacFt.column,"F*Sx*Ft");
798+ printMatrix (tmpmatrixC1.data ,tmpmatrixC1.rows ,tmpmatrixC1.column ," F*Sx" );
799+ printMatrix (tmpmatrixD1.data ,tmpmatrixD1.rows ,tmpmatrixD1.column ," F*Sx*Ft" );
789800 }
790801 /*
791802 * calculate F* for covariance matrix (F*Sx*Ftranspose).F*= (F*Sx)
792- * matrix C will be overridden with new values which is the output
803+ * tmpmatrixC1 will be overridden with new values which is the output
793804 * for the calculation A *x =B
794805 * A = tmpmatrixD
795806 * B = tmpmatrixC
796807 */
797- solveSystemFstar (jacF.rows ,Sx.column ,tmpmatrixD,tmpmatrixC );
808+ solveSystemFstar (jacF.rows ,Sx.column ,tmpmatrixD1. data ,tmpmatrixC1. data );
798809 // printMatrix(tmpmatrixC,jacF.rows,Sx.column,"Sx_F*");
799810 if (ACTIVE_STREAM (LOG_JAC))
800811 {
801- printMatrix (tmpmatrixC ,jacF.rows ,Sx.column ," F*" );
812+ printMatrix (tmpmatrixC1. data ,jacF.rows ,Sx.column ," F*" );
802813 cout << " ***** Completed ****** \n\n " ;
803814 }
804815
805- matrixData tmpFstar = {jacF.rows ,Sx.column ,tmpmatrixC };
816+ matrixData tmpFstar = {jacF.rows ,Sx.column ,tmpmatrixC1. data };
806817 matrixData reconciled_Sx = solveReconciledSx (Sx,jacFt,tmpFstar);
807818 // printMatrix(reconciled_Sx.data,reconciled_Sx.rows,reconciled_Sx.column,"reconciled Sx ===> (Sx - (Sx*Ft*Fstar))");
808819
0 commit comments