Mercurial > forge

changeset 11614:672c998c6620 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
in a mini example / test of pcr usage in the tester. clearing some unused variables declarations.
author michelemartone
date Mon, 08 Apr 2013 16:48:33 +0000
parents e4cb1c2c39c4
children 5d94c1d737da
files main/sparsersb/src/sparsersb.cc main/sparsersb/src/sparsersbtester.m
diffstat 2 files changed, 53 insertions(+), 21 deletions(-) [+]
line wrap: on
line diff
--- a/main/sparsersb/src/sparsersb.cc	Mon Apr 08 11:52:52 2013 +0000
+++ b/main/sparsersb/src/sparsersb.cc	Mon Apr 08 16:48:33 2013 +0000
@@ -310,7 +310,7 @@
 			Array<rsb_coo_idx_t> IA( dim_vector(1,sm.nnz()) );
 			Array<rsb_coo_idx_t> JA( dim_vector(1,sm.nnz()) );
 			rsb_err_t errval=RSB_ERR_NO_ERROR;
-			bool islowtri=true,isupptri=true;
+			/* bool islowtri=sm.is_lower_triangular(),isupptri=sm.is_upper_triangular(); */
 			rsb_flags_t eflags=RSBOI_RF;
 			rsb_type_t typecode=RSB_NUMERICAL_TYPE_DOUBLE;
 			octave_idx_type nrA = sm.rows (), ncA = sm.cols ();
@@ -342,7 +342,7 @@
 			Array<rsb_coo_idx_t> IA( dim_vector(1,sm.nnz()) );
 			Array<rsb_coo_idx_t> JA( dim_vector(1,sm.nnz()) );
 			rsb_err_t errval=RSB_ERR_NO_ERROR;
-			bool islowtri=true,isupptri=true;
+			/* bool islowtri=sm.is_lower_triangular(),isupptri=sm.is_upper_triangular(); */
 			rsb_flags_t eflags=RSBOI_RF;
 			rsb_type_t typecode=RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX;
 #if RSBOI_WANT_SYMMETRY
@@ -356,10 +356,10 @@
 				RSBOI_0_ERROR(RSBOI_0_ALLERRMSG);
 		}
 
-		octave_sparsersb_mtx (const ComplexMatrix &m) : octave_sparse_matrix (RSBIO_DEFAULT_CORE_MATRIX)
+		octave_sparsersb_mtx (const ComplexMatrix &cm) : octave_sparse_matrix (RSBIO_DEFAULT_CORE_MATRIX)
 		{
 			RSBOI_DEBUG_NOTICE(RSBOI_D_EMPTY_MSG);
-			this->alloc_rsb_mtx_from_csc_copy(SparseComplexMatrix(m));
+			this->alloc_rsb_mtx_from_csc_copy(SparseComplexMatrix(cm));
 		}
 
 		octave_sparsersb_mtx (const SparseComplexMatrix &sm, rsb_type_t typecode=RSBOI_TYPECODE) : octave_sparse_matrix (RSBIO_DEFAULT_CORE_MATRIX)
@@ -543,10 +543,11 @@
 					{
 						RSBOI_DEBUG_NOTICE(RSBOI_D_EMPTY_MSG);
 	    					idx_vector i = idx.front() (0).index_vector ();
+#if   defined(RSB_LIBRSB_VER) && (RSB_LIBRSB_VER< 10100)
 						octave_idx_type ii=i(0);
-#if   defined(RSB_LIBRSB_VER) && (RSB_LIBRSB_VER< 10100)
 						RSBOI_ERROR("");
 #elif defined(RSB_LIBRSB_VER) && (RSB_LIBRSB_VER>=10100)
+						octave_idx_type ii=i(0);
 						RSBOI_DEBUG_NOTICE("get_element (%d)\n",ii);
 						if(is_real_type())
 						{
@@ -800,7 +801,6 @@
 				default:
 					panic_impossible ();
 			}
-skipimpl:
 			return retval;
 		}
 
@@ -897,7 +897,7 @@
 			for(nzi=0;nzi<nnzA;++nzi)
 				octave_stdout<<"  ("<<1+IA(nzi)<<", "<<1+JA(nzi)<<") -> "<<((RSBOI_T*)rcm.VA)[nzi]<<"\n";
 			newline(os);
-done:			RSBIO_NULL_STATEMENT_FOR_COMPILER_HAPPINESS
+			RSBIO_NULL_STATEMENT_FOR_COMPILER_HAPPINESS
 		}
 
 	octave_value diag (octave_idx_type k) const
@@ -916,14 +916,12 @@
 			if(this->is_real_type())
 			{
 				Matrix DA(this->rows(),1);
-				//errval=rsb_getdiag (this->mtxAp,(RSBOI_T*)DA.data());
 				errval = rsb_mtx_get_vec(this->mtxAp,(RSBOI_T*)DA.data(),RSB_EXTF_DIAG);
 				retval=(DA);
 			}
 			else
 			{
 				ComplexMatrix DA(this->rows(),1);
-				//errval=rsb_getdiag (this->mtxAp,(void*)DA.data());
 				errval = rsb_mtx_get_vec(this->mtxAp,(RSBOI_T*)DA.data(),RSB_EXTF_DIAG);
 				retval=(DA);
 			}
@@ -937,8 +935,6 @@
 
 	octave_value rsboi_get_scaled_copy_inv(const RSBOI_T alpha)const
 	{
-		rsb_err_t errval=RSB_ERR_NO_ERROR;
-		octave_sparsersb_mtx * m = NULL;
 		RSBOI_T one=1.0;
 		RSBOI_DEBUG_NOTICE(RSBOI_D_EMPTY_MSG);
 		return rsboi_get_scaled_copy(one/alpha);/* FIXME: is this correct ? */
@@ -947,8 +943,6 @@
 #if RSBOI_WANT_DOUBLE_COMPLEX
 	octave_value rsboi_get_scaled_copy_inv(const Complex alpha)const
 	{
-		rsb_err_t errval=RSB_ERR_NO_ERROR;
-		octave_sparsersb_mtx * m = NULL;
 		Complex one=1.0;
 		RSBOI_DEBUG_NOTICE(RSBOI_D_EMPTY_MSG);
 		return rsboi_get_scaled_copy(one/alpha);/* FIXME: is this correct ? */
@@ -958,7 +952,6 @@
 	octave_value rsboi_get_scaled_copy(const RSBOI_T alpha)const
 	{
 		rsb_err_t errval=RSB_ERR_NO_ERROR;
-		octave_sparsersb_mtx * m = NULL;
 		RSBOI_DEBUG_NOTICE(RSBOI_D_EMPTY_MSG);
 		struct rsb_mtx_t*mtxBp=NULL;
 		if(is_real_type())
@@ -974,8 +967,7 @@
 #else
 		{RSBOI_0_ERROR(RSBOI_0_NOCOERRMSG);}
 #endif
-		m = new octave_sparsersb_mtx( mtxBp );
-		return m;
+		return new octave_sparsersb_mtx( mtxBp );
 	}
 
 #if RSBOI_WANT_DOUBLE_COMPLEX
@@ -1009,7 +1001,6 @@
 
 static octave_base_value * default_numeric_conversion_function (const octave_base_value& a)
 {
-	rsb_err_t errval=RSB_ERR_NO_ERROR;
 	RSBOI_DEBUG_NOTICE(RSBOI_D_EMPTY_MSG);
 	CAST_CONV_ARG (const octave_sparsersb_mtx&);
 	RSBOI_WARN(RSBOI_O_MISSIMPERRMSG);
@@ -1020,10 +1011,10 @@
 		return new octave_sparse_complex_matrix (v.sparse_complex_matrix_value());
 }
 
-DEFINE_OCTAVE_ALLOCATOR (octave_sparsersb_mtx);
+DEFINE_OCTAVE_ALLOCATOR (octave_sparsersb_mtx)
 DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparsersb_mtx,
 RSB_OI_TYPEINFO_STRING,
-RSB_OI_TYPEINFO_TYPE);
+RSB_OI_TYPEINFO_TYPE)
 
 DEFCONV (octave_triangular_conv, octave_sparsersb_mtx, matrix)
 {
@@ -1137,9 +1128,13 @@
 	if(RSBOI_SOME_ERROR(errval))
 	{
 		if(errval==RSB_ERR_INVALID_NUMERICAL_DATA)
+		{
 			RSBOI_PERROR(errval);// FIXME: need a specific error message here
+		}
 		else
+		{
 			RSBOI_PERROR(errval);// FIXME: generic case, here
+		}
 		for(octave_idx_type i=0;i<nels;++i)
 			retval(i)=octave_NaN;
 	}
@@ -1158,9 +1153,13 @@
 	if(RSBOI_SOME_ERROR(errval))
 	{
 		if(errval==RSB_ERR_INVALID_NUMERICAL_DATA)
+		{
 			RSBOI_PERROR(errval);// FIXME: need a specific error message here
+		}
 		else
+		{
 			RSBOI_PERROR(errval);// FIXME: generic case, here
+		}
 		for(octave_idx_type i=0;i<nels;++i)
 			retval(i)=octave_NaN;
 	}
@@ -1767,9 +1766,10 @@
 		/* FIXME: undocumented feature */
 		rsb_err_t errval=RSB_ERR_NO_ERROR;
 		const char *mis=args(2).string_value().c_str();
-		rsb_real_t miv=RSBOI_ZERO;/* FIXME: this is extreme danger! */
+		/* rsb_real_t miv=RSBOI_ZERO;*/ /* FIXME: this is extreme danger! */
 		char ss[RSBOI_INFOBUF];
-		if(!osmp || !osmp->mtxAp)goto ret;/* FIXME: error handling missing here */
+		if(!osmp || !osmp->mtxAp)
+			goto ret;/* FIXME: error handling missing here */
 		if(strlen(mis)==0)
 		{
 			mis="RSB_MIF_MATRIX_INFO__TO__CHAR_P";
--- a/main/sparsersb/src/sparsersbtester.m	Mon Apr 08 11:52:52 2013 +0000
+++ b/main/sparsersb/src/sparsersbtester.m	Mon Apr 08 16:48:33 2013 +0000
@@ -200,6 +200,37 @@
 	testmsg(match,"pcg");
 end
 
+function match=testpcrm(OM,XM)
+	# FIXME! This test ignores OM and XM !
+	match=1;
+	A=sparse   ([11,12;21,23]);X=[11;22];B=A*X;X=[0;0];
+	[OX, OFLAG, ORELRES, OITER, ORESVEC]=pcr(A,B);
+	A=sparsersb([11,12;21,23]);X=[11;22];B=A*X;X=[0;0];
+	[XX, XFLAG, XRELRES, XITER, XRESVEC]=pcr(A,B);
+	match&=(OFLAG==XFLAG);# FIXME: a very loose check!
+	match&=(OITER==XITER);# FIXME: a very loose check!
+	#
+	# http://www.gnu.org/software/octave/doc/interpreter/Iterative-Techniques.html#Iterative-Techniques
+	#n = 10;
+	n = 10+size(XM)(1,1)
+	clear A OX XX;
+	A = diag (sparse (1:n));
+	A = A + A'; # we want symmetric matrices
+	b = rand (n, 1);
+	[l, u, p, q] = luinc (A, 1.e-3);
+	[OX, OFLAG, ORELRES, OITER, ORESVEC]= pcr (          A ,b);
+	[XX, XFLAG, XRELRES, XITER, XRESVEC]= pcr (sparsersb(A),b);
+	match&=(sum(OX-XX)==0.0);# FIXME: a very brittle check!
+	#
+	#function y = apply_a (x)
+	#	y = [1:N]' .* x;
+	#endfunction
+	[OX, OFLAG, ORELRES, OITER, ORESVEC]= pcr (          A ,b, 1.e-6, 500, l);
+	[XX, XFLAG, XRELRES, XITER, XRESVEC]= pcr (sparsersb(A),b, 1.e-6, 500, l);
+	match&=(sum(OX-XX)==0.0);# FIXME: a very brittle check!
+	testmsg(match,"pcr");
+end
+
 function match=testmult(OM,XM)
 	match=1;
 	B=ones(rows(OM),1);
@@ -333,6 +364,7 @@
 	match&=testasgn(OM,XM);
 	if nnz(OM)>1
 		match&=testpcgm(OM,XM);
+		match&=testpcrm(OM,XM);
 		match&=testmult(OM,XM);
 		match&=testspsv(OM,XM);
 		match&=testnorm(OM,XM);