Commit b752ac79 authored by Mike Kirby's avatar Mike Kirby

*** empty log message ***


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent c21bdd07
#include "BlockMat.h"
/*
g++ -I../../../include -g -o SCdemo SCdemo.cpp -L../../ -lBlockMat -lblas -llapack -lg2c
*/
using namespace blockmat;
/* declare matrix
| 1 2 - - |
A = | 3 4 - - | using submatrix mat = | 1 2 |
| - - 1 2 | | 3 4 |
| - - 3 4 |
| 1 2 3 - - - - - - |
B = | 4 5 6 - - - - - - | using submatrix mat = | 1 2 3 |
| 1 2 3 1 2 3 1 2 3 | | 4 5 6 |
| 4 5 6 4 5 6 4 5 6 |
| 1 2 3 - - - - - - |
| 4 5 6 - - - - - - |
C = | 7 0 0 - - - - - - | using submatrix mat = | 1 2 3 |
| - - - 1 2 3 - - - | | 4 5 6 |
| - - - 4 5 6 - - - | | 7 0 0 |
| - - - 7 0 0 - - - |
| - - - - - - 1 2 3 |
| - - - - - - 4 5 6 |
| - - - - - - 7 0 0 |
| 1 2 1 2 |
D = | 3 4 3 4 | using submatrix mat = | 1 2 |
| 5 6 5 6 | | 3 4 |
| - - 1 2 | | 5 6 |
| - - 3 4 |
| - - 5 6 |
| - - 1 2 |
| - - 3 4 |
| - - 5 6 |
Calculate DC = A - B*C*D
*/
main(){
BlockMat *A,*B,*C,*D,*T,*SC;
double *mat;
mat = new double [9];
mat[0] = 1; mat[1] = 2; mat[2] = 3;
mat[3] = 4; mat[4] = 5; mat[5] = 6;
mat[6] = 7; mat[7] = 0; mat[8] = 0;
cout << "A: " << endl;
A = new BlockMat(2,2);
A->GenBlock(0,0,2,2,mat);
A->GenBlock(1,1,2,2,mat);
A->PrintBlocks();
cout << endl << "B: " << endl;
B = new BlockMat(2,3);
B->GenBlock(0,0,2,3,mat);
B->GenBlock(1,0,2,3,mat);
B->GenBlock(1,1,2,3,mat);
B->GenBlock(1,2,2,3,mat);
B->PrintBlocks();
cout << endl << "C: " << endl;
C = new BlockMat(3,3);
C->GenBlock(0,0,3,3,mat);
C->GenBlock(1,1,3,3,mat);
C->GenBlock(2,2,3,3,mat);
C->PrintBlocks();
cout << endl << "C^{-1}: " << endl;
C->invert_diag();
C->PrintBlocks();
cout << endl << "D: " << endl;
D = new BlockMat(3,2);
D->GenBlock(0,0,3,2,mat);
D->GenBlock(0,1,3,2,mat);
D->GenBlock(1,1,3,2,mat);
D->GenBlock(2,1,3,2,mat);
D->PrintBlocks();
cout << endl << "SC=A-B*C*D: " << endl;
SC = new BlockMat(2,2);
T = new BlockMat(3,2);
// T->geMxM(RowMajor,RowMajor,1,*C,*D,0);
T->MxM(*C,*D);
SC->sub(*A,SC->MxM(*B,*T));
SC->PrintBlocks();
cout << endl << "SC=A-D^T*C*D: " << endl;
T->MxM(*C,*D);
SC->sub(*A,SC->MtxM(*D,*T));
SC->PrintBlocks();
cout << endl << "SC=A-B*C*B^T: " << endl;
T->MxMt(*C,*B);
//T->geMxM(RowMajor,ColMajor,1,*C,*B,0);
SC->sub(*A,SC->MxM(*B,*T));
SC->PrintBlocks();
delete A;
delete B;
delete C;
delete D;
delete T;
delete SC;
delete[] mat;
return 0;
}
#include "BlockMat.h"
/*
g++ -I../../../include -g -o demo demo.cpp -L../../ -lBlockMat -lblas -lg2c -llapack
*/
using namespace blockmat;
main(){
int i;
BlockMat *A,*B,*C;
double *mat;
/* declare matrix
| 1 2 - - |
A = | 3 4 - - | using submatrix mat = | 1 2 |
| - - 1 2 | | 3 4 |
| - - 3 4 |
*/
mat = new double [4];
mat[0] = 1; mat[1] = 2; mat[2] = 3; mat[3] = 4;
cout << "A: " << endl;
A = new BlockMat(2,2);
A->GenBlock(0,0,2,2,mat);
A->GenBlock(0,1,2,2,mat);
A->GenBlock(1,1,2,2,mat);
A->PrintBlocks();
cout << endl << "B: " << endl;
B = new BlockMat(2,2);
B->GenBlock(0,0,2,2,mat);
B->GenBlock(1,1,2,2,mat);
B->PrintBlocks();
cout << endl << "C=A+B: " << endl;
C = new BlockMat(2,2);
C->add(*A,*B);
C->PrintBlocks();
cout << endl << "C=A-B: " << endl;
C->sub(*A,*B);
C->PrintBlocks();
cout << endl << "C=A*B: " << endl;
C->MxM(*A,*B);
C->PrintBlocks();
double *y = new(double)[4];
double *v = new(double)[4];
Vmath::fill(4,1.0,v,1);
Vmath::zero(4,y,1);
cout << endl << "y = A*v: "<< endl;
C->Mxvpy(v,y);
for(i = 0; i < 4; ++i)
cout << y[i] << " ";
cout << endl;
Vmath::zero(4,y,1);
cout << endl << "y = A^T*v: "<< endl;
//C->Mtxvpy(v,y);
C->geMxv(ColMajor,1,v,1,y);
for(i = 0; i < 4; ++i)
cout << y[i] << " ";
cout << endl;
delete A;
delete B;
delete C;
delete[] y;
delete[] v;
delete[] mat;
return 0;
}
#include <cstdio>
#include <cstdlib>
#include <LibUtilities/NekGraph.h>
int main(int argc, char *argv[])
{
return 0;
}
bin_PROGRAMS = MemoryManagerExample
#GraphExample_SOURCES = GraphExample.cpp
#VecMatExample_SOURCES = VecMatExample.cpp
MemoryManagerExample_SOURCES = MemoryManager.cpp
#GraphExample_CPPFLAGS = -I$(srcdir)/../..
#VecMatExample_CPPFLAGS = -I$(srcdir)/../.. -I$(srcdir)/../../../ThirdParty/met -I$(srcdir)/../../../ThirdParty/met/common
MemoryManagerExample_CPPFLAGS = -I$(srcdir)/../..
MemoryManagerExample_LDADD = -lboost_thread-gcc-mt
This diff is collapsed.
#include <LibUtilities/ThreadSpecificPool.hpp>
#include <LibUtilities/NekMemoryManager.hpp>
#include <boost/thread/thread.hpp>
#include <boost/thread/mutex.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/mpl/size_t.hpp>
#include <iostream>
using namespace std;
using namespace Nektar;
class Disabled
{
public:
Disabled()
{
cout << "Creating a disabled." << endl;
}
~Disabled()
{
cout << "Destroying a disabled." << endl;
}
static const MemoryManager::MemoryPoolEnabler MemoryPoolEnabled = MemoryManager::eDisabled;
};
class Enabled
{
public:
Enabled()
{
cout << "Creating a Enabled." << endl;
}
~Enabled()
{
cout << "Destroying a Enabled." << endl;
}
static const MemoryManager::MemoryPoolEnabler MemoryPoolEnabled = MemoryManager::eEnabled;
};
int main()
{
// Object allocation.
Disabled* t1 = MemoryManager::Allocate<Disabled>();
MemoryManager::Deallocate(t1);
assert(t1 == NULL);
const Disabled* t2 = MemoryManager::Allocate<Disabled>();
MemoryManager::Deallocate(t2);
assert(t2 == NULL);
Enabled* t3 = MemoryManager::Allocate<Enabled>();
MemoryManager::Deallocate(t3);
assert(t3 == NULL);
boost::shared_ptr<Disabled> t4 = MemoryManager::AllocateSharedPtr<Disabled>();
boost::shared_ptr<Enabled> t5 = MemoryManager::AllocateSharedPtr<Enabled>();
// This doesn't compile, it really doesn't make sense to allow it to.
//double* d = MemoryManager::Allocate<double>();
cout << "\nDouble Array." << endl;
double* d = MemoryManager::AllocateArray<10, double>();
MemoryManager::DeallocateArray<10>(d);
cout << "\nDisabled Array" << endl;
Disabled* a1 = MemoryManager::AllocateArray<2, Disabled>();
MemoryManager::DeallocateArray<2>(a1);
cout << "\nEnabled Array." << endl;
Enabled* a2 = MemoryManager::AllocateArray<3, Enabled>();
MemoryManager::DeallocateArray<3>(a2);
boost::shared_array<Disabled> a3 = MemoryManager::AllocateSharedArray<3, Disabled>();
boost::shared_array<Enabled> a4 = MemoryManager::AllocateSharedArray<4, Enabled>();
boost::shared_array<double> a5 = MemoryManager::AllocateSharedArray<5, double>();
boost::shared_array<double> a6 = MemoryManager::AllocateSharedArray<double>(12);
cout << "\nDestroy all shared pointers." << endl;
// boost::shared_array<Class2> a = MemoryManager::AllocateArray<Class2>();
// boost::shared_array<double> d = MemoryManager::AllocateArray<double>();
// Allocate a raw
// boost::thread_group g;
// for(unsigned int i = 0; i < 2; ++i)
// {
// //g.create_thread(testThread);
// g.create_thread(testNakedThread);
// }
//
// g.join_all();
return 0;
}
#include <cstdio>
#include <cstdlib>
#include <LibUtilities/NekLinAlg.h>
const int size = 10;
using namespace LibUtilities::NekLinAlg;
int main(int argc, char *argv[])
{
return 0;
// int i,j;
// double dx = 2.0*M_PI/(double)(size-1);
//
// // Array of doubles
// boost::shared_array<double> dblptr1(new double[size]), dblptr2(new double[size]);
//
// // NekPoint statically allocated
// NekPoint<double,size> pt1,pt2;
//
// // NekPoint dynamically allocated
// NekPoint<double,size> * pt3 = new NekPoint<double,size>,
// *pt4 = new NekPoint<double,size>;
//
// // NekVector statically allocated
// NekVector<double> vec1(size),vec2(size);
//
// // NekVector dynamically allocated
// NekVector<double> * vec3 = new NekVector<double>(size),
// *vec4 = new NekVector<double>(size);
//
// // NekMatrix staticallyy allocated
// NekMatrix<double> mat1(size),mat2(size,size);
//
// // NekMatrix dynamically allocated
// NekMatrix<double> *mat3 = new NekMatrix<double>(size),
// *mat4 = new NekMatrix<double>(size,size);
//
// // Arrays to populate variables
// for(i=0;i<size;i++){
// pt1(1) = 1.0;
// pt2(i) = 1.0 + i*dx;
// (*pt3)(i) = sin(i*dx);
// (*pt4)(i) = cos(i*dx);
// vec1(i) = 1.0;
// vec2(i) = 1.0 + i*dx;
// (*vec3)(i) = sin(i*dx);
// (*vec4)(i) = cos(i*dx);
//
// for(j=0;j<size;j++){
// mat1(i,j) = mat2(i,j) = (*mat3)(i,j) = (*mat4)(i,j) = pow(i*dx,j);
// }
//
// }
//
//
// // Add two points together
// (*pt4) = pt1 + pt2;
//
// // Subtract two points
// (*pt4) = pt1 - pt2;
//
// // Scale a point by a number
// (*pt3) = 2.0*(*pt4);
// (*pt4) = (*pt3)*2.0;
//
// // Add a vector to a point
// (*pt3) = pt1 + vec1;
// (*pt4) = vec2 + pt2;
//
// // Subtract a vector from a point
// pt1 = pt2 - vec2;
//
// // Add two vectors together
// (*vec3) = vec1 + vec2;
//
// // Subtract two points
// (*vec4) = vec1 - vec2;
//
// // Scale a point by a number
// (*vec3) = 2.0*(*vec4);
// (*vec4) = (*vec3)*2.0;
//
// // add mutiple vectors together
// vec1 = vec1 + vec2 + (*vec3) + (*vec4);
//
// // add matrices
// mat1 = mat2 + (*mat3);
//
// // subtract matrices
// mat2 = mat1 - (*mat4);
}
bin_PROGRAMS = Project1D Project2D
Project1D_SOURCES = Project1D.cpp
Project2D_SOURCES = Project2D.cpp
Project1D_CPPFLAGS = -I$(srcdir)/../..
Project2D_CPPFLAGS = -I$(srcdir)/../..
#include <cstdio>
#include <cstdlib>
#include <LocalRegions/LocalRegions.hpp>
#include <LocalRegions/SegExp.h>
#include <SpatialDomains/SegGeom.h>
using namespace Nektar;
using namespace StdRegions;
using namespace LocalRegions;
using namespace SpatialDomains;
using namespace std;
// compile using Builds/Demos/StdRegions -> make DEBUG=1 ProjectS1D
// This routine projects a polynomial or trigonmetric functions which
// has energy in all mdoes of the expansions and report an error.
main(int argc, char *argv[]){
int i,j;
int order, nq;
const double *z,*w;
double *sol,*phys,L2_err;
double x[2];
PointsType Qtype;
BasisType btype;
StdExpansion1D *E;
if(argc != 6){
fprintf(stderr,"Usage: Project1D Type order nq x0 x1 \n");
fprintf(stderr,"Where type is an integer value which "
"dictates the basis as:\n");
fprintf(stderr,"\t Ortho_A = 0\n");
fprintf(stderr,"\t Modified_A = 3\n");
fprintf(stderr,"\t Fourier = 6\n");
fprintf(stderr,"\t Lagrange = 7\n");
fprintf(stderr,"\t Legendre = 8\n");
fprintf(stderr,"\t Chebyshev = 9\n");
fprintf(stderr,"Note type = 1,2,4,5 are for higher dimensional basis\n");
exit(1);
}
btype = (BasisType) atoi(argv[1]);
// Check to see that only 1D Expansions are used
if((btype == eOrtho_B)||(btype == eOrtho_B)||
(btype == eModified_B)||(btype == eModified_C))
{
ErrorUtil::Error(ErrorUtil::efatal,"Project1D",
"This basis is for 2 or 3D expansions");
}
order = atoi(argv[2]);
nq = atoi(argv[3]);
x[0] = atof(argv[4]);
x[1] = atof(argv[5]);
sol = new double [nq];
if(btype != eFourier)
{
Qtype = eLobatto;
}
else
{
Qtype = eFourierEvenSp;
}
//-----------------------------------------------
// Define a segment expansion based on basis definition
const BasisKey b1(btype,order, Qtype, nq,0,0);
VertexComponent *vert1 = new VertexComponent(1,0,x);
VertexComponent *vert2 = new VertexComponent(1,0,x+1);
SegGeom *geom = new SegGeom(vert1,vert2);
geom->SetOwnData();
E = new SegExp(b1,geom);
E->SetGeoFac(E->GenGeoFac());
//-----------------------------------------------
//----------------------------------------------
// Define solution to be projected
E->GetZW(0,z,w);
double *xc = new double [nq];
E->GetCoords(&xc);
if(btype != eFourier)
{
for(i = 0; i < nq; ++i)
{
sol[i] = 0.0;
for(j = 0; j < order; ++j)
{
sol[i] += pow(xc[i],j);
}