Matrix square root with Babylonian algorithm, some updates to ADMFD elasticity example

Examples/ADMFD/elastic.cpp

@@ -455,6 +455,7 @@ int main(int argc,char ** argv)
 		const rMatrix iBtB(viBtB,6,6);
 		const rMatrix Curl(vCurl,9,9);
 		const rMatrix I = rMatrix::Unit(3);
+		const rMatrix I9 = rMatrix::Unit(9);
 		
 		//PrintSV(B);
 		std::cout << "B^T*B" << std::endl;
@@ -492,10 +493,10 @@ int main(int argc,char ** argv)
             ttt = Timer();
             //Assemble gradient matrix W on cells
 #if defined(USE_OMP)
-#pragma omp parallel
+//#pragma omp parallel
 #endif
 			{
-				rMatrix N, R, L, T, K(9,9), C(6,6), W1, W2, W3, U,S,V, w, u, v;
+				rMatrix N, R, L, M, T, K(9,9), C(6,6), W1, W2, W3, U,S,V, w, u, v;
 				rMatrix x(3,1), xf(3,1), n(3,1);
 				double area; //area of the face
 				double volume; //volume of the cell
@@ -518,9 +519,11 @@ int main(int argc,char ** argv)
 					 //K += rMatrix::Unit(9)*1.0e-6*K.FrobeniusNorm();
 					 //PrintSV(K);
 					 N.Resize(3*NF,9); //co-normals
+					 //NQ.Resize(3*NF,9);
 					 //T.Resize(3*NF,9); //transversals
 					 R.Resize(3*NF,9); //directions
 					 L.Resize(3*NF,3*NF);
+					 M.Resize(3*NF,3*NF);
 					 L.Zero();
 					 
 					 //A.Resize(3*NF,3*NF);
@@ -540,6 +543,32 @@ int main(int argc,char ** argv)
 						 
 						 N(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose());
 						 
+						 /*
+						 Cell c2 = cell.Neighbour(faces[k]);
+						 
+						 if( c2.isValid() )
+						 {
+							 KTensor(tag_C[c2],K2);
+							 CTensor(tag_C[c2],C2);
+							 //T1 = I.Kronecker(n.Transpose())*K*I.Kronecker(n);
+							 T2 = I.Kronecker(n.Transpose())*K2*I.Kronecker(n);
+							 //TagRealArray tag_G = m->GetTag("REFERENCE_GRADIENT");
+							 Q1 = I9 + I.Kronecker(n)*T2.Invert()*I.Kronecker(n.Transpose())*(K-K2);
+							 //Q2 = I9 + I.Kronecker(n)*T1.Invert()*I.Kronecker(n.Transpose())*(K2-K);
+							 //NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*(I9 - B*Q*((B*Q).Transpose()*B*Q)*(B*Q).Transpose());
+							 Q1 = Q1.Root();
+							 //NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*(I9 - (Q1.Transpose()*B)*((Q1.Transpose()*B).Transpose()*(Q1.Transpose()*B))*(Q1.Transpose()*B).Transpose());
+							 //NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*Q1*(I9 - B*(B.Transpose()*B)*B.Transpose());
+							 //NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*(I9 - B*(B.Transpose()*B)*B.Transpose())*Q1;
+							 NQ(3*k,3*(k+1),0,9) = area*I.Kronecker(n.Transpose())*Q1;
+							 //NQ(3*k,3*(k+1),0,9).Zero();
+						 }
+						 else
+						 {
+							 NQ(3*k,3*(k+1),0,9).Zero();
+						 }
+						 */
+						 
 						 //std::cout << "I\otimes n^T" << std::endl;
 						 //I.Kronecker(n.Transpose()).Print();
 						 
@@ -642,30 +671,61 @@ int main(int argc,char ** argv)
 					 //W2 = L - (L*R*B)*((L*R*B).Transpose()*R*B).CholeskyInvert()*(L*R*B).Transpose();
 					 //W2 = W1.Trace()*(rMatrix::Unit(3*NF) - (R*B)*((R*B).Transpose()*R*B).Invert()*(R*B).Transpose());
 					 
-					 if( true )
+					 if( false )
 					 {
 						 double alpha = 1;
 						 double beta = alpha;
 						 
 						 //R = R*(rMatrix::Unit(9) + B*iBtB*B.Transpose())*0.5;
-						 K += 1.0e-5*K.Trace()*(rMatrix::Unit(9) - B*iBtB*B.Transpose());
+						 //if( cell.GetElementType() == Element::Tet )
+							//K += 1.0e-3*K.Trace()*(rMatrix::Unit(9) - B*iBtB*B.Transpose());
+						 
+						 
+						 W1 = (N*K+alpha*L*R)*((N*K+alpha*L*R).Transpose()*R).PseudoInvert(1.0e-11)*(N*K+alpha*L*R).Transpose();
+						 
+						 //R += N*(rMatrix::Unit(9) - B*iBtB*B.Transpose())*K.Trace()*2*NF/volume;
+						 
+						 W2 = L - (1+beta)*(L*R)*((L*R).Transpose()*R).PseudoInvert(1.0e-11)*(L*R).Transpose();
 						 
-						 W1 = (N*K+alpha*L*R)*((N*K+alpha*L*R).Transpose()*R).Invert()*(N*K+alpha*L*R).Transpose();
-						 W2 = L - (1+beta)*(L*R)*((L*R).Transpose()*R).PseudoInvert()*(L*R).Transpose();
+						 if( cell.GetElementType() == Element::Tet )
+						 {
+							 //R = R*B*iBtB;
+							 //W2 += (L - (L*R)*((L*R).Transpose()*L*R).PseudoInvert(1.0e-11)*(L*R).Transpose());
+							 N = N*B;
+							 R = R*B*iBtB;//*B.Transpose();
+							 
+							 //W1 = (N*C)*((N*C).Transpose()*R).Invert()*(N*C).Transpose();
+							 W2 += 1.0e-5*(L - (L*R)*((L*R).Transpose()*R).PseudoInvert(1.0e-11)*(L*R).Transpose());
+						 }
+						 
+#pragma omp critical
+						{
+						 std::cout << "W1: "; PrintSV(W1);
+						 std::cout << "W2: "; PrintSV(W2);
+						 std::cout << "S : "; PrintSV(W1+W2);
+						}
 					 }
 					 else
 					 {
-						 //W2 = L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
+						 
 						 
 						 N = N*B;
-						 R = R*B*iBtB;
+						 R = R*B;//*iBtB;
+						 
+						 W1 = (N*C*(B.Transpose()*B))*((N*C*(B.Transpose()*B)).Transpose()*R).Invert()*(N*C*(B.Transpose()*B)).Transpose();
+						 W2 = L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
 						 
-						 W1 = (N*C)*((N*C).Transpose()*R).Invert()*(N*C).Transpose();
 						 //double alpha = 1;
 						 //W1 = (N*K+alpha*L*R)*((N*K+alpha*L*R).Transpose()*R).Invert()*(N*K+alpha*L*R).Transpose() - alpha*(L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
 						 
+#pragma omp critical
+						{
+						 std::cout << "W1: "; PrintSV(W1);
+						 std::cout << "W2: "; PrintSV(W2);
+						 std::cout << "S : "; PrintSV(W1+W2);
+						}
 						 //R = R*B*iBtB;
-						 W2 += L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
+						 //W2 = L - (L*R)*((L*R).Transpose()*R).Invert()*(L*R).Transpose();
 						 //W2 = 2*W1.Trace()/(3*NF) *(rMatrix::Unit(3*NF) - (R)*((R).Transpose()*R).Invert()*(R).Transpose());
 						/*	
 						 #pragma omp critical
@@ -1192,8 +1252,8 @@ int main(int argc,char ** argv)
 				
 				//R.GetJacobian().Save("A.mtx",&Text);
 
-				Solver S(Solver::INNER_MPTILU2);
-                //Solver S(Solver::INNER_MPTILUC);
+				//Solver S(Solver::INNER_MPTILU2);
+                Solver S(Solver::INNER_MPTILUC);
 				//Solver S("superlu");
                 S.SetParameter("relative_tolerance", "1.0e-14");
                 S.SetParameter("absolute_tolerance", "1.0e-12");

Source/Headers/inmost_dense.h

@@ -751,9 +751,25 @@ namespace INMOST
 		/// @param ierr Returns error on fail. If ierr is NULL, then throws an exception.
 		///             If *ierr == -1 on input, then prints out information in case of failure.
 		///             In case of failure *ierr = 1, in case of no failure *ierr = 0.
-		/// @return A pair of pseudo-inverse matrix and boolean. If boolean is true,
-		///         then the matrix was inverted successfully.
+		/// @return A pseudo-inverse of the matrix.
 		Matrix<Var, pool_array_t<Var> > PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0, int * ierr = NULL) const;
+		
+		/// Calcuate A^n, where n is some real value.
+		/// @param n Real value.
+		/// @param ierr Returns error on fail. If ierr is NULL, then throws an exception.
+		///             If *ierr == -1 on input, then prints out information in case of failure.
+		///             In case of failure *ierr = 1, in case of no failure *ierr = 0.
+		///             The error may be caused by error in SVD algorithm.
+		/// @return Matrix in power of n.
+		//Matrix<Var, pool_array_t<Var> > Power(INMOST_DATA_REAL_TYPE n = 1, int * ierr = NULL) const;
+		/// Calculate square root of A matrix by Babylonian method.
+		/// @param iter Number of iterations.
+		/// @param tol  Convergence tolerance.
+		/// @param ierr Returns error on fail. If ierr is NULL, then throws an exception.
+		///             If *ierr == -1 on input, then prints out information in case of failure.
+		///             In case of failure *ierr = 1, in case of no failure *ierr = 0.
+		/// @return Square root of a matrix.
+		Matrix<Var, pool_array_t<Var> > Root(INMOST_DATA_ENUM_TYPE iter = 25, INMOST_DATA_REAL_TYPE tol = 1.0e-7, int * ierr = NULL) const;
 		/// Solves the system of equations of the form A*X=B, with A and B matrices.
 		/// Uses Moore-Penrose pseudo-inverse of the matrix A and calculates X = A^+*B.
 		/// @param B Matrix at the right hand side.
@@ -2571,7 +2587,58 @@ namespace INMOST
 		if( ierr ) *ierr = 0;
 		return ret;
 	}
-	
+	template<typename Var>
+	Matrix<Var, pool_array_t<Var> >
+	AbstractMatrix<Var>::Root(INMOST_DATA_ENUM_TYPE iter, INMOST_DATA_REAL_TYPE tol, int *ierr) const
+	{
+		assert(Rows() == Cols());
+		Matrix<Var, pool_array_t<Var> > ret(Cols(),Cols());
+		Matrix<Var, pool_array_t<Var> > ret0(Cols(),Cols());
+		ret.Zero();
+		ret0.Zero();
+		int k = 0;
+		for(k = 0; k < Cols(); ++k) ret(k,k) = ret0(k,k) = 1;
+		while(k < iter)
+		{
+			ret0 = ret;
+			ret = 0.5*(ret + (*this)*ret.Invert());
+			if( (ret - ret0).FrobeniusNorm() < tol ) return ret; 
+		}
+		if( ierr )
+		{
+			if( *ierr == -1 ) std::cout << "Failed to find square root of matrix by Babylonian method" << std::endl;
+			*ierr = 1;
+			return ret;
+		}
+		return ret;
+	}
+	/*
+	template<typename Var>
+	Matrix<Var, pool_array_t<Var> >
+	AbstractMatrix<Var>::Power(INMOST_DATA_REAL_TYPE n, int * ierr) const
+	{
+		Matrix<Var, pool_array_t<Var> > ret(Cols(),Rows());
+		Matrix<Var, pool_array_t<Var> > L,S,iL;
+		bool success = Eigensolver(L,S,iL);
+		if( !success )
+		{
+			if( ierr )
+			{
+				if( *ierr == -1 ) std::cout << "Failed to compute eigenvalue decomposition of the matrix" << std::endl;
+				*ierr = 1;
+				return ret;
+			}
+			else throw MatrixEigensolverFail;
+		}
+        for(INMOST_DATA_ENUM_TYPE k = 0; k < S.Cols(); ++k) S(k,k) = pow(S(k,k),n);
+        if( n >= 0 )
+			ret = U*S*V.Transpose();
+		else
+			ret = V*S*U.Transpose();
+		if( ierr ) *ierr = 0;
+		return ret;
+	}
+	*/
 	template<typename Var>
 	template<typename typeB>
 	Matrix<typename Promote<Var,typeB>::type, pool_array_t<typename Promote<Var,typeB>::type> >