Creation of a subdirectory mex in the matlab folder (sources + mex files).

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1495 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/mex/sources/gensylv/cc/BlockDiagonal.cpp

+/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/BlockDiagonal.cpp,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
+
+/* Tag $Name:  $ */
+
+#include "BlockDiagonal.h"
+
+#include <stdio.h>
+#include <string.h>
+
+BlockDiagonal::BlockDiagonal(const double* d, int d_size)
+	: QuasiTriangular(d, d_size),
+	  row_len(new int[d_size]), col_len(new int[d_size])
+{
+	for (int i = 0; i < d_size; i++) {
+		row_len[i] = d_size;
+		col_len[i] = 0;
+	}
+}
+
+BlockDiagonal::BlockDiagonal(const QuasiTriangular& t)
+	: QuasiTriangular(t),
+	  row_len(new int[t.numRows()]), col_len(new int[t.numRows()])
+{
+	for (int i = 0; i < t.numRows(); i++) {
+		row_len[i] = t.numRows();
+		col_len[i] = 0;
+	}
+}
+
+BlockDiagonal::BlockDiagonal(int p, const BlockDiagonal& b)
+	: QuasiTriangular(p, b),
+	  row_len(new int[b.numRows()]), col_len(new int[b.numRows()])
+{
+	memcpy(row_len, b.row_len, b.numRows()*sizeof(int));
+	memcpy(col_len, b.col_len, b.numRows()*sizeof(int));
+}
+
+BlockDiagonal::BlockDiagonal(const BlockDiagonal& b)
+	: QuasiTriangular(b),
+	  row_len(new int[b.numRows()]), col_len(new int[b.numRows()])
+{
+	memcpy(row_len, b.row_len, b.numRows()*sizeof(int));
+	memcpy(col_len, b.col_len, b.numRows()*sizeof(int));
+}
+
+/* put zeroes to right upper submatrix whose first column is defined
+ * by 'edge' */
+void BlockDiagonal::setZerosToRU(diag_iter edge)
+{
+	int iedge = (*edge).getIndex();
+	for (int i = 0; i < iedge; i++)
+		for (int j = iedge; j < numCols(); j++)
+			get(i,j) = 0.0;
+}
+
+/* Updates row_len and col_len so that there are zeroes in upper right part, this
+ * |T1 0 |
+ * |0  T2|. The first column of T2 is given by diagonal iterator 'edge'.
+
+ * Note the semantics of row_len and col_len. row_len[i] is distance
+ * of the right-most non-zero element of i-th row from the left, and
+ * col_len[j] is distance of top-most non-zero element of j-th column
+ * to the top. (First element has distance 1).
+ */
+void BlockDiagonal::setZeroBlockEdge(diag_iter edge)
+{
+	setZerosToRU(edge);
+
+	int iedge = (*edge).getIndex();
+	for (diag_iter run = diag_begin(); run != edge; ++run) {
+		int ind = (*run).getIndex();
+		if (row_len[ind] > iedge) {
+			row_len[ind] = iedge;
+			if (!(*run).isReal())
+				row_len[ind+1] = iedge;
+		} 
+	}
+	for (diag_iter run = edge; run != diag_end(); ++run) {
+		int ind = (*run).getIndex();
+		if (col_len[ind] < iedge) {
+			col_len[ind] = iedge;
+			if (!(*run).isReal())
+				col_len[ind+1] = iedge;
+		}
+	}
+}
+
+BlockDiagonal::const_col_iter
+BlockDiagonal::col_begin(const DiagonalBlock& b) const
+{
+	int jbar = b.getIndex();
+	int d_size = diagonal.getSize();
+	return const_col_iter(&getData()[jbar*d_size + col_len[jbar]], d_size,
+						  b.isReal(), col_len[jbar]);
+}
+
+BlockDiagonal::col_iter
+BlockDiagonal::col_begin(const DiagonalBlock& b)
+{
+	int jbar = b.getIndex();
+	int d_size = diagonal.getSize();
+	return col_iter(&getData()[jbar*d_size + col_len[jbar]], d_size,
+					b.isReal(), col_len[jbar]);
+}
+
+
+BlockDiagonal::const_row_iter
+BlockDiagonal::row_end(const DiagonalBlock& b) const
+{
+	int jbar = b.getIndex();
+	int d_size = diagonal.getSize();
+	return const_row_iter(&getData()[d_size*row_len[jbar]+jbar], d_size,
+						  b.isReal(), row_len[jbar]);
+}
+
+BlockDiagonal::row_iter
+BlockDiagonal::row_end(const DiagonalBlock& b)
+{
+	int jbar = b.getIndex();
+	int d_size = diagonal.getSize();
+	return row_iter(&getData()[d_size*row_len[jbar]+jbar], d_size,
+					b.isReal(), row_len[jbar]);
+}
+
+int BlockDiagonal::getNumZeros() const
+{
+	int sum = 0;
+	for (int i = 0; i < diagonal.getSize(); i++) {
+		sum += diagonal.getSize() - row_len[i];
+	}
+	return sum;
+}
+
+QuasiTriangular::const_diag_iter
+BlockDiagonal::findBlockStart(const_diag_iter from) const
+{
+	if (from != diag_end()) {
+		++from;
+		while (from != diag_end() &&
+			   col_len[(*from).getIndex()] != (*from).getIndex())
+			++from;
+	}
+	return from;
+}
+
+int BlockDiagonal::getLargestBlock() const
+{
+	int largest = 0;
+	const_diag_iter start = diag_begin();
+	const_diag_iter end = findBlockStart(start);
+	while (start != diag_end()) {
+		int si = (*start).getIndex();
+		int ei = diagonal.getSize();
+		if (end != diag_end())
+			ei = (*end).getIndex();
+		if (largest < ei-si)
+			largest = ei-si;
+		start = end;
+		end = findBlockStart(start);
+	}
+	return largest;
+}
+
+
+void BlockDiagonal::savePartOfX(int si, int ei, const KronVector& x, Vector& work)
+{
+	for (int i = si; i < ei; i++) {
+		ConstKronVector xi(x, i);
+		Vector target(work, (i-si)*xi.length(), xi.length());
+		target = xi;
+	}
+}
+
+void BlockDiagonal::multKronBlock(const_diag_iter start, const_diag_iter end,
+								  KronVector& x, Vector& work) const
+{
+	int si = (*start).getIndex();
+	int ei = diagonal.getSize();
+	if (end != diag_end())
+		ei = (*end).getIndex();
+	savePartOfX(si, ei, x, work);
+
+	for (const_diag_iter di = start; di != end; ++di) {
+		int jbar = (*di).getIndex();
+		if ((*di).isReal()) {
+			KronVector xi(x, jbar);
+			xi.zeros();
+			Vector wi(work, (jbar-si)*xi.length(), xi.length());
+			xi.add(*((*di).getAlpha()), wi);
+			for (const_row_iter ri = row_begin(*di); ri != row_end(*di); ++ri) {
+				int col = ri.getCol();
+				Vector wj(work, (col-si)*xi.length(), xi.length());
+				xi.add(*ri, wj);
+			}
+		} else {
+			KronVector xi(x, jbar);
+			KronVector xii(x, jbar+1);
+			xi.zeros();
+			xii.zeros();
+			Vector wi(work, (jbar-si)*xi.length(), xi.length());
+			Vector wii(work, (jbar+1-si)*xi.length(), xi.length());
+			xi.add(*((*di).getAlpha()), wi);
+			xi.add((*di).getBeta1(), wii);
+			xii.add((*di).getBeta2(), wi);
+			xii.add(*((*di).getAlpha()), wii);
+			for (const_row_iter ri = row_begin(*di); ri != row_end(*di); ++ri) {
+				int col = ri.getCol();
+				Vector wj(work, (col-si)*xi.length(), xi.length());
+				xi.add(ri.a(), wj);
+				xii.add(ri.b(), wj);
+			}
+		}
+	}
+}
+
+void BlockDiagonal::multKronBlockTrans(const_diag_iter start, const_diag_iter end,
+									   KronVector& x, Vector& work) const
+{
+	int si = (*start).getIndex();
+	int ei = diagonal.getSize();
+	if (end != diag_end())
+		ei = (*end).getIndex();
+	savePartOfX(si, ei, x, work);
+
+	for (const_diag_iter di = start; di != end; ++di) {
+		int jbar = (*di).getIndex();
+		if ((*di).isReal()) {
+			KronVector xi(x, jbar);
+			xi.zeros();
+			Vector wi(work, (jbar-si)*xi.length(), xi.length());
+			xi.add(*((*di).getAlpha()), wi);
+			for (const_col_iter ci = col_begin(*di); ci != col_end(*di); ++ci) {
+				int row = ci.getRow();
+				Vector wj(work, (row-si)*xi.length(), xi.length());
+				xi.add(*ci, wj);
+			}
+		} else {
+			KronVector xi(x, jbar);
+			KronVector xii(x, jbar+1);
+			xi.zeros();
+			xii.zeros();
+			Vector wi(work, (jbar-si)*xi.length(), xi.length());
+			Vector wii(work, (jbar+1-si)*xi.length(), xi.length());
+			xi.add(*((*di).getAlpha()), wi);
+			xi.add((*di).getBeta2(), wii);
+			xii.add((*di).getBeta1(), wi);
+			xii.add(*((*di).getAlpha()), wii);
+			for (const_col_iter ci = col_begin(*di); ci != col_end(*di); ++ci) {
+				int row = ci.getRow();
+				Vector wj(work, (row-si)*xi.length(), xi.length());
+				xi.add(ci.a(), wj);
+				xii.add(ci.b(), wj);
+			}
+		}
+	}
+}
+
+void BlockDiagonal::multKron(KronVector& x) const
+{
+	int largest = getLargestBlock();
+	Vector work(largest*x.getN()*power(x.getM(),x.getDepth()-1));
+	const_diag_iter start = diag_begin();
+	const_diag_iter end = findBlockStart(start);
+	while (start != diag_end()) {
+		multKronBlock(start, end, x, work);
+		start = end;
+		end = findBlockStart(start);
+	}
+}
+
+
+void BlockDiagonal::multKronTrans(KronVector& x) const
+{
+	int largest = getLargestBlock();
+	Vector work(largest*x.getN()*power(x.getM(),x.getDepth()-1));
+	const_diag_iter start = diag_begin();
+	const_diag_iter end = findBlockStart(start);
+	while (start != diag_end()) {
+		multKronBlockTrans(start, end, x, work);
+		start = end;
+		end = findBlockStart(start);
+	}
+}
+
+void BlockDiagonal::printInfo() const
+{
+	printf("Block sizes:");
+	int num_blocks = 0;
+	const_diag_iter start = diag_begin();
+	const_diag_iter end = findBlockStart(start);
+	while (start != diag_end()) {
+		int si = (*start).getIndex();
+		int ei = diagonal.getSize();
+		if (end != diag_end())
+			ei = (*end).getIndex();
+		printf(" %d", ei-si);
+		num_blocks++;
+		start = end;
+		end = findBlockStart(start);
+	}
+	printf("\nNum blocks: %d\n", num_blocks);
+	printf("There are %d zeros out of %d\n",
+		   getNumZeros(), getNumOffdiagonal());
+}
+
+int BlockDiagonal::getNumBlocks() const
+{
+	int num_blocks = 0;
+	const_diag_iter start = diag_begin();
+	const_diag_iter end = findBlockStart(start);
+	while (start != diag_end()) {
+		num_blocks++;
+		start = end;
+		end = findBlockStart(start);
+	}
+	return num_blocks;
+}
+
+// Local Variables:
+// mode:C++
+// End:

matlab/mex/sources/gensylv/cc/BlockDiagonal.h

+/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/BlockDiagonal.h,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
+
+/* Tag $Name:  $ */
+
+#ifndef BLOCK_DIAGONAL_H
+#define BLOCK_DIAGONAL_H
+
+#include "QuasiTriangular.h"
+
+
+class BlockDiagonal : public QuasiTriangular {
+	int* const row_len;
+	int* const col_len;
+public:
+	BlockDiagonal(const double* d, int d_size);
+	BlockDiagonal(int p, const BlockDiagonal& b);
+	BlockDiagonal(const BlockDiagonal& b);
+	BlockDiagonal(const QuasiTriangular& t);
+	const BlockDiagonal& operator=(const QuasiTriangular& t)
+		{GeneralMatrix::operator=(t); return *this;}
+	const BlockDiagonal& operator=(const BlockDiagonal& b);
+	~BlockDiagonal() {delete [] row_len; delete [] col_len;}
+	void setZeroBlockEdge(diag_iter edge);
+	int getNumZeros() const;
+	int getNumBlocks() const;
+	int getLargestBlock() const;
+	void printInfo() const;
+
+	void multKron(KronVector& x) const;
+	void multKronTrans(KronVector& x) const;
+
+	const_col_iter col_begin(const DiagonalBlock& b) const;
+	col_iter col_begin(const DiagonalBlock& b);
+	const_row_iter row_end(const DiagonalBlock& b) const;
+	row_iter row_end(const DiagonalBlock& b);
+	QuasiTriangular* clone() const
+		{return new BlockDiagonal(*this);}
+private:
+	void setZerosToRU(diag_iter edge);
+	const_diag_iter findBlockStart(const_diag_iter from) const;
+	static void savePartOfX(int si, int ei, const KronVector& x, Vector& work);
+	void multKronBlock(const_diag_iter start, const_diag_iter end,
+					   KronVector& x, Vector& work) const;
+	void multKronBlockTrans(const_diag_iter start, const_diag_iter end,
+							KronVector& x, Vector& work) const;
+};
+
+#endif /* BLOCK_DIAGONAL_H */
+
+
+// Local Variables:
+// mode:C++
+// End:

matlab/mex/sources/gensylv/cc/GeneralMatrix.cpp

+/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralMatrix.cpp,v 1.4 2004/11/24 20:41:59 kamenik Exp $ */
+
+/* Tag $Name:  $ */
+
+
+#include "SylvException.h"
+#include "GeneralMatrix.h"
+
+#include "cppblas.h"
+#include "cpplapack.h"
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <cmath>
+#include <limits>
+
+int GeneralMatrix::md_length = 32;
+
+GeneralMatrix::GeneralMatrix(const GeneralMatrix& m)
+	: data(m.rows*m.cols), rows(m.rows), cols(m.cols), ld(m.rows)
+{
+	copy(m);
+}
+
+GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix& m)
+	: data(m.rows*m.cols), rows(m.rows), cols(m.cols), ld(m.rows)
+{
+	copy(m);
+}
+
+GeneralMatrix::GeneralMatrix(const GeneralMatrix& m, const char* dummy)
+	: data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols)
+{
+	for (int i = 0; i < m.rows; i++)
+		for (int j = 0; j < m.cols; j++)
+			get(j,i) = m.get(i,j);
+}
+
+GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix& m, const char* dummy)
+	: data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols)
+{
+	for (int i = 0; i < m.rows; i++)
+		for (int j = 0; j < m.cols; j++)
+			get(j,i) = m.get(i,j);
+}
+
+
+GeneralMatrix::GeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols)
+	: data(nrows*ncols), rows(nrows), cols(ncols), ld(nrows)
+{
+	copy(m, i, j);
+}
+
+GeneralMatrix::GeneralMatrix(GeneralMatrix& m, int i, int j, int nrows, int ncols)
+	: data(m.base()+m.ld*j+i, m.ld*(ncols-1)+nrows), rows(nrows), cols(ncols), ld(m.ld)
+{}
+
+GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b)
+	: data(a.rows*b.cols), rows(a.rows), cols(b.cols), ld(a.rows)
+{
+	gemm("N", a, "N", b, 1.0, 0.0);
+}
+
+GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b, const char* dum)
+	: data(a.rows*b.rows), rows(a.rows), cols(b.rows), ld(a.rows)
+{
+	gemm("N", a, "T", b, 1.0, 0.0);
+}
+
+GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b)
+	: data(a.cols*b.cols), rows(a.cols), cols(b.cols), ld(a.cols)
+{
+	gemm("T", a, "N", b, 1.0, 0.0);
+}
+
+GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum1,
+							 const GeneralMatrix& b, const char* dum2)
+	: data(a.cols*b.rows), rows(a.cols), cols(b.rows), ld(a.cols)
+{
+	gemm("T", a, "T", b, 1.0, 0.0);
+}
+
+
+
+GeneralMatrix::~GeneralMatrix()
+{
+}
+
+
+
+void GeneralMatrix::place(const ConstGeneralMatrix& m, int i, int j)
+{
+	if (i + m.numRows() > numRows() ||
+		j + m.numCols() > numCols())
+		throw SYLV_MES_EXCEPTION("Bad submatrix placement, matrix dimensions exceeded.");
+
+	GeneralMatrix tmpsub(*this, i, j, m.numRows(), m.numCols());
+	tmpsub.copy(m);
+}
+
+void GeneralMatrix::mult(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b)
+{
+	gemm("N", a, "N", b, 1.0, 0.0);
+}
+
+void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
+							   double mult)
+{
+	gemm("N", a, "N", b, mult, 1.0);
+}
+
+void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
+							   const char* dum, double mult)
+{
+	gemm("N", a, "T", b, mult, 1.0);
+}
+
+void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const char* dum,
+							   const ConstGeneralMatrix& b, double mult)
+{
+	gemm("T", a, "N", b, mult, 1.0);
+}
+
+void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const char* dum1,
+							   const ConstGeneralMatrix& b, const char* dum2, double mult)
+{
+	gemm("T", a, "T", b, mult, 1.0);
+}
+
+void GeneralMatrix::multRight(const ConstGeneralMatrix& m)
+{
+	gemm_partial_right("N", m, 1.0, 0.0);
+}
+
+void GeneralMatrix::multLeft(const ConstGeneralMatrix& m)
+{
+	gemm_partial_left("N", m, 1.0, 0.0);
+}
+
+void GeneralMatrix::multRightTrans(const ConstGeneralMatrix& m)
+{
+	gemm_partial_right("T", m, 1.0, 0.0);
+}
+
+void GeneralMatrix::multLeftTrans(const ConstGeneralMatrix& m)
+{
+	gemm_partial_left("T", m, 1.0, 0.0);
+}
+
+// here we must be careful for ld
+void GeneralMatrix::zeros()
+{
+	if (ld == rows)
+		data.zeros();
+	else {
+		for (int i = 0; i < rows; i++) 
+			for (int j = 0; j < cols; j++)
+				get(i,j) = 0;
+	}
+}
+
+void GeneralMatrix::unit()
+{
+	for (int i = 0; i < rows; i++)
+		for (int j = 0; j < cols; j++)
+			if (i == j)
+				get(i,j) = 1.0;
+			else
+				get(i,j) = 0.0;
+}
+
+void GeneralMatrix::nans()
+{
+	for (int i = 0; i < rows; i++) 
+		for (int j = 0; j < cols; j++)
+			get(i,j) = std::numeric_limits<double>::quiet_NaN();
+}
+
+void GeneralMatrix::infs()
+{
+	for (int i = 0; i < rows; i++) 
+		for (int j = 0; j < cols; j++)
+			get(i,j) = std::numeric_limits<double>::infinity();
+}
+
+
+// here we must be careful for ld
+void GeneralMatrix::mult(double a)
+{
+	if (ld == rows)
+		data.mult(a);
+	else {