Commit 334063c7 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

sync (cmake download superlu, block of matrix, fix modifications)

parent 30570c97
._*
.svn
Tests/solver_test001/matrices
SyncToy*
......
......@@ -35,6 +35,7 @@ option(USE_SOLVER_MONDRIAAN "Use Mondriaan for matrix reordering" OFF)
option(USE_SOLVER_PETSC "Use PETSc solvers" OFF)
option(USE_SOLVER_TRILINOS "Use Trilinos solvers" OFF)
option(USE_SOLVER_SUPERLU "Use SuperLU solver" OFF)
option(USE_SOLVER_SUPERLU_DOWNLOAD "Attempts to download SuperLU solver if library not found" OFF)
#option(USE_AUTODIFF_ASMJIT "Use AsmJit for automatic differentiation" OFF)
#option(USE_AUTODIFF_EXPRESSION_TEMPLATES "Use c++ expression templates for automatic differentiation" OFF)
......@@ -245,8 +246,56 @@ endif()
if(USE_SOLVER_SUPERLU)
find_package(SUPERLU)
if(NOT SUPERLU_FOUND)
set(USE_SOLVER_SUPERLU OFF CACHE BOOL "Use SuperLU solver" FORCE)
#set(USE_SOLVER_SUPERLU OFF CACHE BOOL "Use SuperLU solver" FORCE)
MESSAGE("SUPERLU NOT FOUND")
if( USE_SOLVER_SUPERLU_DOWNLOAD )
include(ExternalProject)
message("-- Download SUPERLU")
#string (REPLACE ";" "|" PASS_MPI "${MPI_INCLUDE_PATH}")
#message("-- Passing MPI include to ParMETIS: '${PASS_MPI}'")
ExternalProject_Add(SUPERLU
GIT_REPOSITORY "https://github.com/xiaoyeli/superlu.git"
GIT_TAG "master"
UPDATE_DISCONNECTED 1
#${EXTERNAL_NOUPDATE}
PREFIX "${LIB_DOWNLOAD_PATH}"
#this command is responsible for update of repository from git
#UPDATE_COMMAND ""
LIST_SEPARATOR | # Use the alternate list separator
#for the one from ibaned
#CMAKE_ARGS "-DMPI_INCLUDE_PATH:PATH='${PASS_MPI}'"
CMAKE_ARGS "-DCMAKE_INSTALL_PREFIX=${LIB_DOWNLOAD_PATH} "
"-DCMAKE_INSTALL_INCLUDEDIR=${LIB_DOWNLOAD_PATH}/include/superlu "
"-DCMAKE_BUILD_TYPE=Release "
"-DMPI_C_COMPILER=${MPI_C_COMPILER} "
"-DMPI_CXX_COMPILER=${MPI_CXX_COMPILER} "
"-DCMAKE_POSITION_INDEPENDENT_CODE=ON "
"-DCMAKE_C_COMPILER=${CMAKE_C_COMPILER} "
"-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER} "
)
find_package(BLAS)
if (NOT BLAS_FOUND)
if( WIN32 )
set(BLAS_LIBRARIES "${LIB_DOWNLOAD_PATH}/src/SUPERLU-build/CBLAS/blas.lib" CACHE FILEPATH "BLAS_LIBRARIES" FORCE)
else( WIN32 )
set(BLAS_LIBRARIES "${LIB_DOWNLOAD_PATH}/src/SUPERLU-build/CBLAS/libblas.a" CACHE FILEPATH "BLAS_LIBRARIES" FORCE)
endif( WIN32 )
endif()
set(SUPERLU_INCLUDES "${LIB_DOWNLOAD_PATH}/include" CACHE PATH "SUPERLU_INCLUDE_DIR is set" FORCE)
if( WIN32 )
set(SUPERLU_LIBRARIES "${LIB_DOWNLOAD_PATH}/lib/superlu.lib;${BLAS_LIBRARIES}" CACHE PATH "SUPERLU_LIBRARIES is set" FORCE)
else( WIN32 )
set(SUPERLU_LIBRARIES "${LIB_DOWNLOAD_PATH}/lib/libsuperlu.a;${BLAS_LIBRARIES}" CACHE PATH "SUPERLU_LIBRARIES is set" FORCE)
endif( WIN32 )
message("-- SUPERLU_INCLUDE_DIR: ${SUPERLU_INCLUDES}")
message("-- SUPERLU_LIBRARIES: ${SUPERLU_LIBRARIES}")
include_directories(${SUPERLU_INCLUDES})
else(USE_SOLVER_SUPERLU_DOWNLOAD)
set(USE_SOLVER_SUPERLU OFF CACHE BOOL "Use SuperLU solver" FORCE)
endif(USE_SOLVER_SUPERLU_DOWNLOAD)
else()
MESSAGE("\nFound SuperLU! Here are the details: ")
MESSAGE("INCLUDES: ${SUPERLU_INCLUDES}")
......
......@@ -386,6 +386,8 @@ int main(int argc,char ** argv)
TagRealArray tag_W; // Gradient matrix
TagRealArray tag_Ws; // Matrix to reconstruct stress from fluxes
TagRealArray tag_FLUX; // Flux (for error)
TagRealArray tag_f_coefs, tag_i_coefs, tag_f_rhs, tag_i_rhs;
TagReferenceArray tag_f_elems, tag_i_elems;
TagIntegerArray tag_Row;
if( m->GetProcessorsNumber() > 1 ) //skip for one processor job
......@@ -491,6 +493,12 @@ int main(int argc,char ** argv)
tag_BC = m->GetTag("BOUNDARY_CONDITION");
//initialize unknowns at boundary
}
tag_i_elems = m->CreateTag("i_elems",DATA_REFERENCE,FACE,NONE);
tag_f_elems = m->CreateTag("f_elems",DATA_REFERENCE,FACE,NONE);
tag_i_coefs = m->CreateTag("i_coefs",DATA_REAL,FACE,NONE);
tag_f_coefs = m->CreateTag("f_coefs",DATA_REAL,FACE,NONE);
tag_i_rhs = m->CreateTag("i_rhs",DATA_REAL,FACE,NONE,3);
tag_f_rhs = m->CreateTag("f_rhs",DATA_REAL,FACE,NONE,3);
tag_W = m->CreateTag("W",DATA_REAL,CELL,NONE);
tag_Ws = m->CreateTag("Ws",DATA_REAL,CELL,NONE);
......@@ -613,9 +621,8 @@ int main(int argc,char ** argv)
} //end of loop over cells
/*
rMatrix C1(3,3),C2(3,3), iC12sum, C1sum(3,3), C2sum(3,3), CC(3,3);
raMatrix W[2];
rMatrix C1(3,3),C2(3,3), iC12sum, C1sum(3,3), C2sum(3,3), CC(3,3), CCI(3,3), CCF(3,3), B_D(3,3), B_N(3,3), B_R(3,1);
rMatrix W[2];
#if defined(USE_OMP)
#pragma omp for
#endif
......@@ -634,14 +641,14 @@ int main(int argc,char ** argv)
faces[q] = cell->getFaces(); //obtain faces of the cell
W[q] = raMatrixMake(cell->RealArray(tag_W).data(),3*faces[q].size(),3*faces[q].size());
}
Storage::reference_array f_elems = f->ReferenceArray(elems);
Storage::reference_array i_elems = f->ReferenceArray(intrp_elems);
Storage::real_array f_coefs = f->RealArray(trans);
Storage::real_array i_coefs = f->RealArray(intrp_trans);
real & f_rhs = f->Real(rhs);
real & i_rhs = f->Real(intrp_rhs);
Storage::reference_array f_elems = f->ReferenceArray(tag_f_elems);
Storage::reference_array i_elems = f->ReferenceArray(tag_i_elems);
Storage::real_array f_coefs = f->RealArray(tag_f_coefs);
Storage::real_array i_coefs = f->RealArray(tag_i_coefs);
tag_f_rhs(f,3,1).Zero();
tag_i_rhs(f,3,1).Zero();
real aF = 1;
f_rhs = i_rhs = 0;
//f_rhs = i_rhs = 0;
real multi = 1, multf = -aF, multfbnd = -aF;
if( cells.size() == 2 ) //internal face
{
......@@ -651,77 +658,93 @@ int main(int argc,char ** argv)
C1sum.Zero();
C2sum.Zero();
//add back cell faces and compute sum
for(int q = 0; q < (int)faces[0].size(); ++q) if( fabs(W[0](rows[0],q)) > 1.0e-12 )
for(int q = 0; q < (int)faces[0].size(); ++q)
{
CC = W[0](rows[0]*3,rows[0]*3+3,q*3+3);
C1sum += CC;
CC = iC12sum*CC;
if( faces[0][q] != f )
CC = W[0](rows[0]*3,rows[0]*3+3,q*3,q*3+3);
if( CC.FrobeniusNorm() > 1.0e-12 )
{
f_elems.push_back(faces[0][q]);
//f_coefs.push_back(-W[0](rows[0],q)*W2/(W1+W2)*multf);
f_coefs.push_back(-CC*C2*multf);
i_elems.push_back(faces[0][q]);
i_coefs.push_back(-CC*multi);
C1sum += CC;
CC = iC12sum*CC;
if( faces[0][q] != f )
{
CCI = -CC*multi;
CCF = -CC*C2*multf;
i_coefs.insert(i_coefs.end(),CCI.data(),CCI.data()+9);
f_coefs.insert(f_coefs.end(),CCF.data(),CCF.data()+9);
i_elems.push_back(faces[0][q]);
f_elems.push_back(faces[0][q]);
}
}
}
//add front cell faces and compute sum
for(int q = 0; q < (int)faces[1].size(); ++q) if( fabs(W[1](rows[1],q)) > 1.0e-12 )
for(int q = 0; q < (int)faces[1].size(); ++q)
{
CC = W[1](rows[1]*3,rows[1]*3+3,q*3+3);
C2sum += CC;
CC = iC12sum*CC;
if( faces[1][q] != f )
CC = W[1](rows[1]*3,rows[1]*3+3,q*3,q*3+3);
if( CC.FrobeniusNorm() > 1.0e-12 )
{
f_elems.push_back(faces[1][q]);
f_coefs.push_back(CC*C1*multf);
i_elems.push_back(faces[1][q]);
i_coefs.push_back(-CC*multi);
C2sum += CC;
CC = iC12sum*CC;
if( faces[1][q] != f )
{
CCI = -CC*multi;
CCF = CC*C1*multf;
i_coefs.insert(i_coefs.end(),CCI.data(),CCI.data()+9);
f_coefs.insert(f_coefs.end(),CCF.data(),CCF.data()+9);
i_elems.push_back(faces[1][q]);
f_elems.push_back(faces[1][q]);
}
}
}
//already have back cell and front cell
f_coefs[0] = iC12sum*C1sum*C2*multf; //BackCell
raMatrixMake(f_coefs.data()+0,3,3) = iC12sum*C1sum*C2*multf; //BackCell
assert(f_elems[0] == f->BackCell());
f_coefs[1] = -iC12sum*C1*C2sum*multf; //FrontCell
raMatrixMake(f_coefs.data()+9,3,3) = -iC12sum*C1*C2sum*multf; //FrontCell
assert(f_elems[1] == f->FrontCell());
i_coefs[0] = iC12sum*C1sum*multi;
i_coefs[1] = iC12sum*C2sum*multi;
//i_coefs[0] = iC12sum*C1sum*multi;
//i_coefs[1] = iC12sum*C2sum*multi;
}
else //Boundary face
{
real alpha = 0, beta = 1, gamma = 0; //alpha p + beta K \nabla p \cdot n = gamma
if( bnd_conds.isValid() && f->HaveData(bnd_conds) )
B_D.Zero();
B_N = I;
B_R.Zero();
if( tag_BC.isValid() && f.HaveData(tag_BC) )
{
real_array bnd = f->RealArray(bnd_conds);
alpha = bnd[0];
beta = bnd[1];
gamma = bnd[2];
f.UnitNormal(n.data());
GetBC(f.RealArray(tag_BC),n,B_D,B_N,B_R);
}
real Wsum = 0.0, Wc = W[0](rows[0],rows[0])*aF;
f_rhs = 0;
//f_rhs = gamma*Wc/(alpha+beta*Wc)*aF;
i_rhs = gamma/(alpha+beta*Wc);
C1 = W[0](rows[0]*3,rows[0]*3+3,rows[0]*3,rows[0]*3+3)*aF;
C1sum.Zero();
//TODO
//iC12sum = (alpha*I+beta*C1).Invert();
tag_f_rhs(f,3,1).Zero();
tag_i_rhs(f,3,1) = iC12sum*B_R;
for(int q = 0; q < (int)faces[0].size(); ++q) if( fabs(W[0](rows[0],q)) > 1.0e-12 )
{
f_elems.push_back(faces[0][q]);
f_coefs.push_back(-W[0](rows[0],q)*multfbnd);
if( faces[0][q] != f )
CC = W[0](rows[0]*3,rows[0]*3+3,q*3,q*3+3);
if( CC.FrobeniusNorm() > 1.0e-12 )
{
//TODO
//f_elems.push_back(faces[0][q]);
//f_coefs.push_back(W[0](rows[0],q)*alpha/(alpha+beta*Wc)*aF);
i_elems.push_back(faces[0][q]);
i_coefs.push_back(-W[0](rows[0],q)*beta/(alpha+beta*Wc));
//f_coefs.push_back(-CC*multfbnd);
if( faces[0][q] != f )
{
//i_elems.push_back(faces[0][q]);
//i_coefs.push_back(-beta*iC12sum*CC);
}
C1sum += CC;
}
Wsum += W[0](rows[0],q);
}
//f_coefs[0] = -alpha/(alpha+beta*Wc)*Wsum*aF;
f_coefs[0] = Wsum*multfbnd;
assert(f_elems[0] == f->BackCell());
i_coefs[0] = beta*Wsum/(alpha+beta*Wc);
//TODO
//f_coefs[0] = C1sum*multfbnd;
//assert(f_elems[0] == f->BackCell());
//i_coefs[0] = beta*C1sum*iC12sum;
}
} //end of loop over faces
*/
}
std::cout << "Construct W matrix: " << Timer() - ttt << std::endl;
......
......@@ -488,7 +488,7 @@ namespace INMOST
void svg_draw_faces_nc(std::ostream & file, std::vector<face2gl> & set, bool drawedges, double modelview[16], double projection[16], int viewport[4], int highlight)
{
file << "<g stroke=\"none\" fill=\"green\" fill-opacity=\"0.1\">" << std::endl;
file << "<g stroke=\"none\" fill=\"green\" fill-opacity=\"0.1\" stroke-opacity=\"0.1\">" << std::endl;
for (INMOST_DATA_ENUM_TYPE q = 0; q < set.size(); q++) set[q].svg_draw(file, drawedges, modelview, projection, viewport);
file << "</g>" << std::endl;
if (highlight != -1)
......
#ifndef _INC_GLUT_H
#define _INC_GLUT_H
#define MAC_WORKAROUND
#if defined (__APPLE__) || defined(MACOSX)
#include <GLUT/glut.h>
#endif
......
......@@ -260,7 +260,10 @@ void reshape(int w, int h)
width = w;
height = h;
set_matrix3d();
//latest mac os issue
//#if !(defined(MAC_WORKAROUND) && (defined (__APPLE__) || defined(MACOSX)))
glViewport(0, 0, w, h);
//#endf //MAC_WORKAROUND
}
......@@ -271,7 +274,7 @@ double mymy = 0;
void myclickmotion(int nmx, int nmy) // Mouse
{
double lmx = 2.*(nmx/(double)width - 0.5),lmy = 2.*(0.5 - nmy/(double)height), dmx = lmx-mymx, dmy = lmy - mymy;
if( actionstate == 1 ) //middle button
if( actionstate == 1 )//|| (glutGetModifiers() & GLUT_ACTIVE_CTRL) ) //middle button
{
double shiftmod[3] = {0,0,0};
shiftmod[0] += dmx*zoom*std::max(std::max( sright-sleft, stop-sbottom ), sfar-snear);
......@@ -299,7 +302,7 @@ void myclickmotion(int nmx, int nmy) // Mouse
mymx = lmx;
mymy = lmy;
}
else if( actionstate == 2 ) //right button
else if( actionstate == 2 )// || (glutGetModifiers() & GLUT_ACTIVE_SHIFT) ) //right button
{
if( planecontrol )
{
......@@ -2082,7 +2085,7 @@ void svg_draw(std::ostream & file)
}
else
{
std::vector<face2gl> sorted_clip_boundary;
if( oclipper || bclipper )
{
/*
......@@ -2095,7 +2098,7 @@ void svg_draw(std::ostream & file)
svg_draw_faces(file,:drawedges && ::drawedges != 2,temp_boundary,modelview,projection,viewport);
*/
std::vector<face2gl> sorted_clip_boundary(clip_boundary);
sorted_clip_boundary.insert(sorted_clip_boundary.begin(),clip_boundary.begin(),clip_boundary.end());
for(INMOST_DATA_ENUM_TYPE q = 0; q < sorted_clip_boundary.size() ; q++)
sorted_clip_boundary[q].compute_dist(campos);
face2gl::radix_sort_dist(sorted_clip_boundary);
......@@ -2127,15 +2130,42 @@ void svg_draw(std::ostream & file)
for (int k = 0; k < segments.size(); ++k)
segments[k].SVGDraw(file, modelview, projection, viewport);
file << "</g>" << std::endl;
if( boundary )
if( (oclipper || bclipper) && boundary )
{
file << "<g fill=\"black\" fill-opacity=\"0.25\">" << std::endl;
face2gl::radix_sort_dist(all_boundary);
unsigned i = 0, j = 0;
while( i < sorted_clip_boundary.size() || j < all_boundary.size() )
{
if( j == all_boundary.size() || (i != sorted_clip_boundary.size() && all_boundary[j] < sorted_clip_boundary[i]) )
{
sorted_clip_boundary[i].svg_draw_colour(file, ::drawedges && ::drawedges != 2, modelview, projection, viewport);
i++;
}
else
{
file << "<g stroke=\"none\" fill=\"green\" fill-opacity=\"0.1\" stroke-opacity=\"0.3\">" << std::endl;
all_boundary[j].svg_draw(file, ::drawedges && ::drawedges != 2, modelview, projection, viewport);
file << "</g>" << std::endl;
j++;
}
}
}
else if( oclipper || bclipper )
{
if (!(drawedges == 2 || drawedges == 3))
svg_draw_faces(file, sorted_clip_boundary, ::drawedges && ::drawedges != 2, modelview, projection, viewport);
}
else if( boundary )
{
//file << "<g fill=\"black\" fill-opacity=\"0.25\">" << std::endl;
//if (drawedges && drawedges != 2)svg_draw_edges(file,all_boundary,modelview,projection,viewport);
if (!(drawedges == 2 || drawedges == 3))
svg_draw_faces_nc(file, all_boundary, ::drawedges && ::drawedges != 2, modelview, projection, viewport);
file << "</g>" << std::endl;
//file << "</g>" << std::endl;
}
}
......
#include "inc_glut.h"
#include "tga.h"
#include <stdio.h>
extern void draw_screen(); //global openl drawing routine
......@@ -15,7 +15,9 @@ void screenshot(int tiles)
char * pixelbuffer = new char[width*height * 3 + oldwidth*oldheight * 3];
char * tempbuffer = pixelbuffer + width*height * 3;
//int WH[2];
//glGetIntegerv(GL_VIEWPORT,WH);
//printf("W %d H %d\n",WH[0],WH[1]);
for (int i = 0; i < tiles; ++i)
......@@ -57,5 +59,7 @@ void screenshot(int tiles)
delete[] pixelbuffer;
width = oldwidth;
height = oldheight;
//#if !(defined(MAC_WORKAROUND) && (defined (__APPLE__) || defined(MACOSX)))
glViewport(0, 0, width, height);
//#endf //MAC_WORKAROUND
}
......@@ -111,7 +111,7 @@ static bool BoundingEllipse(Element e, double eps, int iters, rMatrix & Q, rMatr
{
for (int k = 0; k < p.Rows(); ++k) D(k, k) = p(k, 0);
M = A1.Transpose()*(A1*D*A1.Transpose()).Invert()*A1; // m by m
M = A1.Transpose()*(A1*D*A1.Transpose()).Solve(A1); // m by m
//std::cout << "matrix M:" << std::endl;
//M.Print();
kp = -1.0e20;
......@@ -161,7 +161,7 @@ static bool BoundingEllipse(Element e, double eps, int iters, rMatrix & Q, rMatr
}
for (int k = 0; k < p.Rows(); ++k) D(k, k) = p(k, 0);
Q = (A*(D - p*p.Transpose())*A.Transpose()).Invert() / static_cast<double>(d);
Q = (A*(D - p*p.Transpose())*A.Transpose()).PseudoInvert() / static_cast<double>(d);
c = A*p;
//check
if (d == 2)
......@@ -405,6 +405,18 @@ int main(int argc, char ** argv)
std::cout << "collapse edges: " << nedges << std::endl;
*/
/*
for(int k = 0; k < m.EdgeLastLocalID(); ++k) if( m.isValidEdge(k) )
{
Edge e = m.EdgeByLocalID(k);
if( e.GetMarker(collapse) && e.Boundary() )
{
std::cout << "Unmark boundary edge " << e.LocalID() << std::endl;
e.RemMarker(collapse);
}
}
*/
/*
std::vector<HandleType> edges;
for(int k = 0; k < m.EdgeLastLocalID(); ++k) if( m.isValidEdge(k) )
{
......@@ -413,16 +425,32 @@ int main(int argc, char ** argv)
edges.push_back(e.GetHandle());
}
std::sort(edges.begin(),edges.end(),Mesh::MeasureComparator(&m));
*/
int ncollapsed = 0;
m.Save("collapsed.pmf");
//scanf("%*c");
for(unsigned k = 0; k < edges.size(); ++k)
bool found_edge;
do
{
if( !m.isValidEdge(GetHandleID(edges[k])) ) continue; //it may be already deleted
Node n = Edge::CollapseEdge(Edge(&m,edges[k]),0);
cosm[n] = 1;
std::cout << n.LocalID() << std::endl;
ncollapsed++;
found_edge = false;
Edge e = InvalidEdge();
for(Mesh::iteratorEdge it = m.BeginEdge(); it != m.EndEdge(); ++it)
if( it->GetMarker(collapse) )
{
if( e == InvalidEdge() || e.Length() > it->Length() )
{
e = it->self();
found_edge = true;
}
}
if( e.isValid() )
{
std::cout << "collapse edge " << e.LocalID() << std::endl;
Node n = Edge::CollapseEdge(e,0);
cosm[n] = 1;
std::cout << "new node " << n.LocalID() << std::endl;
ncollapsed++;
}
m.Save("collapsed.pmf");
//scanf("%*c");
......@@ -433,9 +461,10 @@ int main(int argc, char ** argv)
throw -1;
}
}
while( found_edge );
std::cout << "Collapsed: " << ncollapsed << std::endl;
std::cout << "Time to fix tiny:" << Timer() - tt << std::endl;
std::cout << "Cells: " << m.NumberOfCells() << std::endl;
std::cout << "Faces: " << m.NumberOfFaces() << std::endl;
......
This diff is collapsed.
......@@ -54,7 +54,7 @@
// output xml files for debugging of parallel algorithms
// search for style.xsl within examples for comfortable
// view of generated xml files
//#define USE_PARALLEL_WRITE_TIME
#define USE_PARALLEL_WRITE_TIME
// this will revert Mesh::PrepareReceiveInner to always
// use MPI point to point functionality disregarding problem type
......@@ -286,6 +286,12 @@ namespace INMOST
template<typename Var>
class ConstSubMatrix;
template<typename Var>
class BlockOfMatrix;
template<typename Var>
class ConstBlockOfMatrix;
template<typename Var, typename Storage = array<Var> >
class Matrix;
......
......@@ -538,7 +538,7 @@ namespace INMOST
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type, pool_array_t<typename Promote<Var,typeB>::type> >
CrossProduct(const AbstractMatrix<typeB> & other) const;
/// Transformation matrix from current vector to provided vector.
/// Transformation matrix from current vector to provided vector using shortest arc rotation.
/// @param other Vector to transform to.
/// @return A sqaure (rotation) matrix that transforms current vector into right hand side vector.
/// \warning Works only for 3x1 vectors.
......@@ -824,6 +824,18 @@ namespace INMOST
SubMatrix<Var> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
ConstSubMatrix<Var> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col) const;
/// Define matrix as a part of a matrix of larger size with in-place manipulation of elements.
/// Let A = {a_ij}, i in [0,n), j in [0,m) be original matrix.
/// Then the method returns B = {a_ij}, if i in [offset_row,offset_row+n),
/// and j in [offset_col,offset_col+m) and B = {0} otherwise.
/// @param nrows Number of rows in larger matrix.
/// @param ncols Number of columns in larger matrix.
/// @param offset_row Offset for row number.
/// @param offset_col Offset for column number.
/// @return Submatrix of the original matrix.
BlockOfMatrix<Var> BlockOf(enumerator nrows, enumerator ncols, enumerator offset_row, enumerator offset_col);
ConstBlockOfMatrix<Var> BlockOf(enumerator nrows, enumerator ncols, enumerator offset_row, enumerator offset_col) const;
};
......@@ -2001,6 +2013,197 @@ namespace INMOST
}
};
/// This class allows to address a matrix as a block of an empty matrix of larger size.
template<typename Var>
class BlockOfMatrix : public AbstractMatrix<Var>
{
public:
using AbstractMatrix<Var>::operator();
typedef unsigned enumerator; //< Integer type for indexes.
private:
AbstractMatrix<Var> * M;
enumerator nrows; //< Number of rows in larger matrix.
enumerator ncols; //< Number of colums in larger matrix.
enumerator orow; //< Row offset in larger matrix.
enumerator ocol; //< Column offset in larger matrix.
public:
/// Number of rows in submatrix.
/// @return Number of rows.
__INLINE enumerator Rows() const {return nrows;}
/// Number of columns in submatrix.
/// @return Number of columns.
__INLINE enumerator Cols() const {return ncols;}
/// Create submatrix for a matrix.
/// @param rM Reference to the matrix that stores elements.
/// @param num_rows Number of rows in the larger matrix.
/// @param num_cols Number of columns in the larger matrix.
/// @param first_row Offset for row index in the larger matrix.
/// @param first_column Offset for column index in the larger matrix.
BlockOfMatrix(AbstractMatrix<Var> & rM, enumerator num_rows, enumerator num_cols, enumerator offset_row, enumerator offset_col) : M(&rM), nrows(num_rows), ncols(num_cols), orow(offset_row), ocol(offset_col)
{}
BlockOfMatrix(const BlockOfMatrix & b) : M(b.M), nrows(b.nrows), ncols(b.ncols), orow(b.orow), ocol(b.ocol) {}
/// Assign matrix of another type to the block of matrix.
/// \warning Only part of the matrix related to non-empty block is copied, other entries are ignored.
/// @param other Another matrix of different type.
/// @return Reference to current submatrix.
template<typename typeB>
BlockOfMatrix & operator =(AbstractMatrix<typeB> const & other)
{
assert( Cols() == other.Cols() );
assert( Rows() == other.Rows() );
for(enumerator i = orow; i < orow+M->Rows(); ++i)
for(enumerator j = ocol; j < ocol+M->Cols(); ++j)
assign((*M)(i-orow,j-ocol),other(i,j));
return *this;
}
/// Assign matrix of another type to the block of matrix.
/// \warning Only part of the matrix related to non-empty block is copied, other entries are ignored.
/// @return Reference to current submatrix.
template<typename typeB>
BlockOfMatrix & operator =(BlockOfMatrix<typeB> const & other)
{
assert( Cols() == other.Cols() );
assert( Rows() == other.Rows() );
for(enumerator i = orow; i < orow+M->Rows(); ++i)
for(enumerator j = ocol; j < ocol+M->Cols(); ++j)