Commit 09c09f3b by Kirill Terekhov

add some files

parent f823b8e7
......@@ -2,12 +2,8 @@ projects/
AsmJit/
compile32_vs2013/
compile64_vs2013/
development/
inmost_autodiff.h
autodiff.cpp
TODO.txt
examples/
examples/Data/
examples/Matrices/
examples/Grids/
tests/
\ No newline at end of file
tests/
//FOR 2D
// Node -> Vertex
// Edge -> Vertex
// Face -> Line, Curve
// Cell -> Tri, Quad, Polygon, MultiLine
//FOR 3D
// Node -> Vertex
// Edge -> Line, Curve
// Face -> Tri, Quad, Polygon, MultiLine
// Cell -> Tet, Hex, Prism, Pyramid, Polyhedron, MultiPolygon
//TODO:
1D???
//serial
//1.0 read/write formats
// 1.1 read gmv
// 1.2 read different kinds of vtk besides unstructured
// 1.3 read/write xml vtk formats
// 1.4 pmf: print error when markers cannot be represented in current type and data was lost
// 1.5 pmf: save hide_marker, new_marker values
// 1.6 pmf: save topology checking information
//2. enchance iterators
// iterators should return pointers or references?
//3. geometry.cpp
// 3.0 decide how FixNormalOrientation should update geometric data
// 3.1 replace calculation of normale by hyperplane equation for any dimension??
// 3.2 correct Inside function to check weather point inside Segment, Polygon
// 3.3 dynamically balanced obtree or octree or kd-tree to query elements by position
//4. (done, need check) add some low-level mesh modification procedures to elements
// 5.1 add some high-level mesh modification procedures: split by plane,
// 5.2 csg: difference, union, intersection
// 5.3 calculate volume for concave elements, choose elements by minimum volume, not minimum number of elements in incident_matrix class (modify.cpp)
//5. reduce size used by sparse data
//6. new ElementSet class, derived from Element class
// 6.1 support ierarhy
//7. algorithm that checks topology on element creation
// 7.1 complete unimplemented tests for star-shaped objects, concavity
//parallel:
//1. algorithm in EndModification or ResolveModification to keep the mesh in parallel consistent state after element creation/deletion
// 1.1 restore contiguous global id
//2. test algorithm that checks for halo in ExchangeMarked, if performance increases
//3. decide how to exchange ElementSets between processors
//4. exchange data of type DATA_REFERENCE, by GlobalID and ElementType
//5. mesh_file.cpp
// 5.1 on parallel load, when there are more processors then parts in mesh - duplicate mesh on other processors??
// 5.2 when mesh is duplicated over many processors use local clusterization algorithm in ResolveShared
// shared parallel
//1. avoid markers in nbAdjElements, getAdjElements, getNodes, getEdges, getFaces, getCells, BridgeAdjacency
// partitioner:
//1. implement K-means clustering
//solver:
//0. workaround for overlapping vertices and matrices
//1. read TODO on top of solver_*.hpp
//autodiff:
//1. read TODO on top of INMOST_autodiff.h
//TODO:
// Recheck computation of derivatives!!!!
// Detect similar parts of the tree in Automatizator
// Restructure tree expressions
// Intorduce multivariate tables
// Generate opencl code
// Make so that stencil may be represented by tags, by set or by callback_function
//RegisterTable and table implementation should account for boundary conditions beyond data
//SOLVER todo:
//TODO:
// how to calculate diagonal perturbation?
// how to calculate schur complement faster
//
// done! implement crout-ILU from
// Documents/Read/solver/crout/crout.pdf to ILUC_preconditioner
// done! implement condition estimation for |L^{-1}| and |U^{-1}| and adaptive tau for ILUC_preconditioner from
// Documents\Read\solver\Read_now\bollhofer_siam.pdf
// done! try to make ILUC2_preconditioner - second order ILU from ILUC_preconditioner
// done! implement diagonal pivoting for ILUC - maybe need to update diagonal at every step
// goto references [7],[10]-(data structure!) from
// Documents\Read\solver\crout\crout.pdf
// return dropped values to diagonal if control vector is provided from
// Documents\Read\solver\stabilization\milut.pdf
// try to apply dropping while forming linked list,should correctly calculate condition estimates
// Calculate schur complement faster:
// Documents\Read\solver\sparse_matmul\sparse.pdf
// in ILUC_preconditioner, replace matrix structures by CSR, estimate number of nonzeros in rows/cols
// before filling, if necessery (how?)
This diff is collapsed. Click to expand it.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
#pragma once
#ifndef __SOLVER_RESCALE__
#define __SOLVER_RESCALE__
//TODO:
//test implementation
//implement RSS scaling from Documents/Read/solver/scaling/2007KaporinRJNAMM.pdf
#include "inmost_solver.h"
#include "solver_prototypes.hpp"
using namespace INMOST;
class RAS_rescale : public Method //Row-Alternating Scaling for 2-norm
{
Solver::Matrix * Alink;
Solver::Vector * xlink;
Solver::Vector * blink;
Solver::OrderInfo * info;
Solver::Vector DL,DR;
INMOST_DATA_ENUM_TYPE iters;
INMOST_DATA_REAL_TYPE subst, eps;
bool init;
public:
RAS_rescale(Solver::Matrix & A, Solver::Vector & x, Solver::Vector & b, Solver::OrderInfo & Minfo)
:Alink(&A),xlink(&x),blink(&b),info(&Minfo)
{
iters = 5;
eps = 1.0e-54;
subst = 1.0;
}
void Copy(const Method * other)
{
const RAS_rescale * b = dynamic_cast<const RAS_rescale *>(other);
assert(b != NULL);
Alink = b->Alink;
info = b->info;
DL = b->DL;
DR = b->DR;
iters = b->iters;
eps = b->eps;
subst = b->subst;
}
RAS_rescale(const RAS_rescale & other)
{
Copy(&other);
}
RAS_rescale & operator =(RAS_rescale const & other)
{
Copy(&other);
return *this;
}
INMOST_DATA_ENUM_TYPE & EnumParameter(std::string name)
{
if( name == "iters" )
return iters;
throw -1;
}
INMOST_DATA_REAL_TYPE & RealParameter(std::string name)
{
if( name == "substitute" )
return subst;
else if( name == "epsilon" )
return eps;
throw -1;
}
bool Initialize()
{
if (isInitialized()) Finalize();
INMOST_DATA_ENUM_TYPE mobeg,moend, vlocbeg, vlocend, vbeg, vend,k, end, iter,r;
info->GetOverlapRegion(info->GetRank(),mobeg,moend);
info->GetLocalRegion(info->GetRank(),vlocbeg,vlocend);
info->GetVectorRegion(vbeg,vend);
DL.SetInterval(mobeg,moend);
info->PrepareVector(DR);
for(k = mobeg; k < moend; k++)
{
for(Solver::Row::iterator q = (*Alink)[k].Begin(); q != (*Alink)[k].End(); ++q)
DL[k] += q->second*q->second;
if( DL[k] < eps ) DL[k] = 1.0/subst; else DL[k] = 1.0/DL[k];
}
for(iter = 0; iter < iters; iter++)
{
std::fill(DR.Begin(),DR.End(),0.0);
for(k = vlocbeg; k < vlocend; k++)
for(Solver::Row::iterator q = (*Alink)[k].Begin(); q != (*Alink)[k].End(); ++q) DR[q->first] += DL[k]*q->second*q->second;
info->Accumulate(DR);
info->Update(DR);
for(k = vbeg; k < vend; k++)
if( DR[k] < eps ) DR[k] = 1.0/subst; else DR[k] = 1.0/DR[k];
std::fill(DL.Begin(),DL.End(),0.0);
for(k = mobeg; k < moend; k++)
for(Solver::Row::iterator q = (*Alink)[k].Begin(); q != (*Alink)[k].End(); ++q) DL[k] += DR[q->first]*q->second*q->second;
for(k = mobeg; k < moend; k++)
if( DL[k] < eps ) DL[k] = 1.0/subst; else DL[k] = 1.0/DL[k];
}
for(k = mobeg; k < moend; k++) DL[k] = sqrt(DL[k]);
for(k = vbeg; k < vend; k++) DR[k] = sqrt(DR[k]);
//Rescale matrix
for(k = mobeg; k < moend; k++)
{
for(Solver::Row::iterator q = (*Alink)[k].Begin(); q != A[k].End(); ++q)
q->second = DL[k] * q->second * DR[q->first];
}
for(k = mobeg; k < moend; k++) (*blink)[k] *= DL[k];
info->Update(*blink);
for(k = vbeg; k < vend; k++) (*xlink)[k] /= DR[k];
init = true;
return true;
}
bool isInitialized() { return init; }
bool ReplaceMAT(Solver::Matrix & A) { if (isInitialized()) Finalize(); Alink = &A; return true; }
bool ReplaceRHS(Solver::Vector & b)
{
INMOST_DATA_ENUM_TYPE mobeg,moend,k;
info->GetOverlapRegion(info->GetRank(),mobeg,moend);
//restore old rhs
for(k = mobeg; k < moend; k++) (*blink)[k] /= DL[k];
//select new rhs
blink = &b;
//rescale new rhs
for(k = mobeg; k < moend; k++) (*blink)[k] *= DL[k];
info->Update(*blink);
return true;
}
bool ReplaceSOL(Solver::Vector & x)
{
INMOST_DATA_ENUM_TYPE vbeg,vend,k;
info->GetVectorRegion(vbeg,vend);
//restore old sol
for(k = vbeg; k < vend; k++) (*xlink)[k] *= DR[k];
//select new sol
xlink = &x;
//rescale new sol
for(k = vbeg; k < vend; k++) (*xlink)[k] /= DR[k];
return true;
}
bool Solve(Solver::Vector & x, Solver::Vector & out)
{
// A x = b
// DL A DR DR^{-1} x = DL b
//rescale vector x by DR
INMOST_DATA_ENUM_TYPE vbeg,vend;
info->GetVectorRegion(vbeg,vend);
for(INMOST_DATA_ENUM_TYPE k = vbeg; k != vend; k++) out[k] = x[k] * DR[k];
return true;
}
bool Finalize()
{
if (!isFinalized())
{
INMOST_DATA_ENUM_TYPE mobeg, moend, k, end, r, vbeg, vend;
info->GetOverlapRegion(info->GetRank(), mobeg, moend);
info->GetVectorRegion(vbeg, vend);
for (k = mobeg; k < moend; k++)
{
Solver::Row & Ak = (*Alink)[k];
end = Ak.Size();
for (r = 0; r < end; r++)
Ak.GetValue(r) = Ak.GetValue(r) / DL[k] / DR[Ak.GetIndex(r)];
}
for (k = mobeg; k < moend; k++) (*blink)[k] /= DL[k];
info->Update(*blink);
for (k = vbeg; k < vend; k++) (*xlink)[k] *= DR[k];
init = false;
}
return true;
}
bool isFinalized() { return !init; }
~RAS_rescale()
{
if (!isFinalized()) Finalize();
}
Method * Duplicate() { return new RAS_rescale(*this); }
};
#endif
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment