DataVectorOps.h
Go to the documentation of this file.
/*****************************************************************************
*
* Copyright (c) 2003-2017 by The University of Queensland
* http://www.uq.edu.au
*
* Primary Business: Queensland, Australia
* Licensed under the Apache License, version 2.0
* http://www.apache.org/licenses/LICENSE-2.0
*
* Development until 2012 by Earth Systems Science Computational Center (ESSCC)
* Development 2012-2013 by School of Earth Sciences
* Development from 2014 by Centre for Geoscience Computing (GeoComp)
*
*****************************************************************************/
#ifndef __ESCRIPT_DATAMATHS_H__
#define __ESCRIPT_DATAMATHS_H__
#include "DataAbstract.h"
#include "DataException.h"
#include "ArrayOps.h"
#include "LapackInverseHelper.h"
#include "DataTagged.h"
#include <complex>
namespace escript
{
ESCRIPT_DLL_API
void
matMult(const DataTypes::RealVectorType& left,
const DataTypes::ShapeType& leftShape,
DataTypes::RealVectorType::size_type leftOffset,
const DataTypes::RealVectorType& right,
const DataTypes::ShapeType& rightShape,
DataTypes::RealVectorType::size_type rightOffset,
DataTypes::RealVectorType& result,
const DataTypes::ShapeType& resultShape);
// Hmmmm why is there no offset for the result??
ESCRIPT_DLL_API
DataTypes::ShapeType
determineResultShape(const DataTypes::ShapeType& left,
const DataTypes::ShapeType& right);
template<typename VEC>
inline
void
symmetric(const VEC& in,
const DataTypes::ShapeType& inShape,
typename VEC::size_type inOffset,
VEC& ev,
const DataTypes::ShapeType& evShape,
typename VEC::size_type evOffset)
{
if (DataTypes::getRank(inShape) == 2) {
int i0, i1;
int s0=inShape[0];
int s1=inShape[1];
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] + in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
}
}
}
else if (DataTypes::getRank(inShape) == 4) {
int i0, i1, i2, i3;
int s0=inShape[0];
int s1=inShape[1];
int s2=inShape[2];
int s3=inShape[3];
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)] + in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)]) / 2.0;
}
}
}
}
}
}
template<typename VEC>
inline
void
antisymmetric(const VEC& in,
const DataTypes::ShapeType& inShape,
typename VEC::size_type inOffset,
VEC& ev,
const DataTypes::ShapeType& evShape,
typename VEC::size_type evOffset)
{
if (DataTypes::getRank(inShape) == 2) {
int i0, i1;
int s0=inShape[0];
int s1=inShape[1];
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] - in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
}
}
}
else if (DataTypes::getRank(inShape) == 4) {
int i0, i1, i2, i3;
int s0=inShape[0];
int s1=inShape[1];
int s2=inShape[2];
int s3=inShape[3];
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)] - in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)]) / 2.0;
}
}
}
}
}
}
void
hermitian(const DataTypes::CplxVectorType& in,
const DataTypes::ShapeType& inShape,
DataTypes::CplxVectorType::size_type inOffset,
DataTypes::CplxVectorType& ev,
const DataTypes::ShapeType& evShape,
DataTypes::CplxVectorType::size_type evOffset);
void
antihermitian(const DataTypes::CplxVectorType& in,
const DataTypes::ShapeType& inShape,
typename DataTypes::CplxVectorType::size_type inOffset,
DataTypes::CplxVectorType& ev,
const DataTypes::ShapeType& evShape,
typename DataTypes::CplxVectorType::size_type evOffset);
template <class VEC>
inline
void
trace(const VEC& in,
const DataTypes::ShapeType& inShape,
typename VEC::size_type inOffset,
VEC& ev,
const DataTypes::ShapeType& evShape,
typename VEC::size_type evOffset,
int axis_offset)
{
for (int j=0;j<DataTypes::noValues(evShape);++j)
{
ev[evOffset+j]=0;
}
if (DataTypes::getRank(inShape) == 2) {
int s0=inShape[0]; // Python wrapper limits to square matrix
int i;
for (i=0; i<s0; i++) {
ev[evOffset/*+DataTypes::getRelIndex(evShape)*/] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
}
}
else if (DataTypes::getRank(inShape) == 3) {
if (axis_offset==0) {
int s0=inShape[0];
int s2=inShape[2];
int i0, i2;
for (i0=0; i0<s0; i0++) {
for (i2=0; i2<s2; i2++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
}
}
}
else if (axis_offset==1) {
int s0=inShape[0];
int s1=inShape[1];
int i0, i1;
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
}
}
}
}
else if (DataTypes::getRank(inShape) == 4) {
if (axis_offset==0) {
int s0=inShape[0];
int s2=inShape[2];
int s3=inShape[3];
int i0, i2, i3;
for (i0=0; i0<s0; i0++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
}
}
}
}
else if (axis_offset==1) {
int s0=inShape[0];
int s1=inShape[1];
int s3=inShape[3];
int i0, i1, i3;
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
}
}
}
}
else if (axis_offset==2) {
int s0=inShape[0];
int s1=inShape[1];
int s2=inShape[2];
int i0, i1, i2;
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
}
}
}
}
}
}
ESCRIPT_DLL_API
template <class VEC>
inline
void
transpose(const VEC& in,
const DataTypes::ShapeType& inShape,
typename VEC::size_type inOffset,
VEC& ev,
const DataTypes::ShapeType& evShape,
typename VEC::size_type evOffset,
int axis_offset)
{
int inRank=DataTypes::getRank(inShape);
if ( inRank== 4) {
int s0=evShape[0];
int s1=evShape[1];
int s2=evShape[2];
int s3=evShape[3];
int i0, i1, i2, i3;
if (axis_offset==1) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
}
}
}
}
}
else if (axis_offset==2) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
}
}
}
}
}
else if (axis_offset==3) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
}
}
}
}
}
else {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
}
}
}
}
}
}
else if (inRank == 3) {
int s0=evShape[0];
int s1=evShape[1];
int s2=evShape[2];
int i0, i1, i2;
if (axis_offset==1) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
}
}
}
}
else if (axis_offset==2) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
}
}
}
}
else {
// Copy the matrix unchanged
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
}
}
}
}
}
else if (inRank == 2) {
int s0=evShape[0];
int s1=evShape[1];
int i0, i1;
if (axis_offset==1) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
}
}
}
else {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
}
}
}
}
else if (inRank == 1) {
int s0=evShape[0];
int i0;
for (i0=0; i0<s0; i0++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
}
}
else if (inRank == 0) {
ev[evOffset/*+DataTypes::getRelIndex(evShape,)*/] = in[inOffset/*+DataTypes::getRelIndex(inShape,)*/];
}
else {
throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
}
}
ESCRIPT_DLL_API
template <class VEC>
inline
void
swapaxes(const VEC& in,
const DataTypes::ShapeType& inShape,
typename VEC::size_type inOffset,
VEC& ev,
const DataTypes::ShapeType& evShape,
typename VEC::size_type evOffset,
int axis0,
int axis1)
{
int inRank=DataTypes::getRank(inShape);
if (inRank == 4) {
int s0=evShape[0];
int s1=evShape[1];
int s2=evShape[2];
int s3=evShape[3];
int i0, i1, i2, i3;
if (axis0==0) {
if (axis1==1) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
}
}
}
}
} else if (axis1==2) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
}
}
}
}
} else if (axis1==3) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
}
}
}
}
}
} else if (axis0==1) {
if (axis1==2) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
}
}
}
}
} else if (axis1==3) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
}
}
}
}
}
} else if (axis0==2) {
if (axis1==3) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
for (i3=0; i3<s3; i3++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
}
}
}
}
}
}
} else if ( inRank == 3) {
int s0=evShape[0];
int s1=evShape[1];
int s2=evShape[2];
int i0, i1, i2;
if (axis0==0) {
if (axis1==1) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
}
}
}
} else if (axis1==2) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
}
}
}
}
} else if (axis0==1) {
if (axis1==2) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
for (i2=0; i2<s2; i2++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
}
}
}
}
}
} else if ( inRank == 2) {
int s0=evShape[0];
int s1=evShape[1];
int i0, i1;
if (axis0==0) {
if (axis1==1) {
for (i0=0; i0<s0; i0++) {
for (i1=0; i1<s1; i1++) {
ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
}
}
}
}
} else {
throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
}
}
ESCRIPT_DLL_API
inline
void
eigenvalues(const DataTypes::RealVectorType& in,
const DataTypes::ShapeType& inShape,
typename DataTypes::RealVectorType::size_type inOffset,
DataTypes::RealVectorType& ev,
const DataTypes::ShapeType& evShape,
typename DataTypes::RealVectorType::size_type evOffset)
{
typename DataTypes::RealVectorType::ElementType in00,in10,in20,in01,in11,in21,in02,in12,in22;
typename DataTypes::RealVectorType::ElementType ev0,ev1,ev2;
int s=inShape[0];
if (s==1) {
in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
eigenvalues1(in00,&ev0);
ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
} else if (s==2) {
in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
} else if (s==3) {
in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
&ev0,&ev1,&ev2);
ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
}
}
inline
void
eigenvalues(const DataTypes::CplxVectorType& in,
const DataTypes::ShapeType& inShape,
typename DataTypes::CplxVectorType::size_type inOffset,
DataTypes::CplxVectorType& ev,
const DataTypes::ShapeType& evShape,
typename DataTypes::CplxVectorType::size_type evOffset)
{
typename DataTypes::CplxVectorType::ElementType in00,in10,in20,in01,in11,in21,in02,in12,in22;
typename DataTypes::CplxVectorType::ElementType ev0,ev1,ev2;
int s=inShape[0];
if (s==1) {
in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
eigenvalues1(in00,&ev0);
ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
} else if (s==2) {
in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
} else if (s==3) {
// this doesn't work yet
// in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
// in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
// in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
// in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
// in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
// in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
// in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
// in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
// in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
// eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
// &ev0,&ev1,&ev2);
// ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
// ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
// ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
}
}
ESCRIPT_DLL_API
inline
void
eigenvalues_and_eigenvectors(const DataTypes::RealVectorType& in, const DataTypes::ShapeType& inShape,
DataTypes::RealVectorType::size_type inOffset,
DataTypes::RealVectorType& ev, const DataTypes::ShapeType& evShape,
DataTypes::RealVectorType::size_type evOffset,
DataTypes::RealVectorType& V, const DataTypes::ShapeType& VShape,
DataTypes::RealVectorType::size_type VOffset,
const double tol=1.e-13)
{
double in00,in10,in20,in01,in11,in21,in02,in12,in22;
double V00,V10,V20,V01,V11,V21,V02,V12,V22;
double ev0,ev1,ev2;
int s=inShape[0];
if (s==1) {
in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
} else if (s==2) {
in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
&ev0,&ev1,&V00,&V10,&V01,&V11,tol);
ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
} else if (s==3) {
in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
&ev0,&ev1,&ev2,
&V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
}
}
template <class VEC>
inline
bool
checkOffset(const VEC& data,
const DataTypes::ShapeType& shape,
typename VEC::size_type offset)
{
return (data.size() >= (offset+DataTypes::noValues(shape)));
}
template <class ResVEC, class LVEC, class RSCALAR>
void
binaryOpVectorRightScalar(ResVEC& res, // where result is to be stored
typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
const typename ResVEC::size_type sampleSize, // number of values in each sample
const LVEC& left, // LHS of calculation
typename LVEC::size_type leftOffset, // where to start reading LHS values
const RSCALAR* right, // RHS of the calculation
const bool rightreset, // true if RHS is providing a single sample of 1 value only
escript::ES_optype operation, // operation to perform
bool singleleftsample) // set to false for normal operation
{
size_t substep=(rightreset?0:1);
switch (operation)
{
case ADD:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
const RSCALAR* rpos=right+(rightreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=left[leftbase+j]+*rpos;
}
}
}
break;
case POW:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
const RSCALAR* rpos=right+(rightreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],*rpos);
}
}
}
break;
case SUB:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
const RSCALAR* rpos=right+(rightreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=left[leftbase+j]-*rpos;
}
}
}
break;
case MUL:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
const RSCALAR* rpos=right+(rightreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=left[leftbase+j] * *rpos;
}
}
}
break;
case DIV:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
const RSCALAR* rpos=right+(rightreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=left[leftbase+j]/ *rpos;
}
}
}
break;
default:
throw DataException("Unsupported binary operation");
}
}
template<>
void
binaryOpVectorRightScalar(DataTypes::RealVectorType& res, // where result is to be stored
typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
const DataTypes::RealVectorType& left, // LHS of calculation
typename DataTypes::RealVectorType::size_type leftOffset, // where to start reading LHS values
const DataTypes::real_t* right, // RHS of the calculation
const bool rightreset, // true if RHS is providing a single sample of 1 value only
escript::ES_optype operation, // operation to perform
bool singleleftsample);
template <class ResVEC, class LSCALAR, class RVEC>
void
binaryOpVectorLeftScalar(ResVEC& res, // where result is to be stored
typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
const typename ResVEC::size_type sampleSize, // number of values in each sample
const LSCALAR* left, // LHS of calculation
const bool leftreset, // true if LHS is providing a single sample of 1 value only
const RVEC& right, // RHS of the calculation
typename RVEC::size_type rightOffset, // where to start reading RHS values
escript::ES_optype operation, // operation to perform
bool singlerightsample) // right consists of a single sample
{
size_t substep=(leftreset?0:1);
switch (operation)
{
case ADD:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
const LSCALAR* lpos=left+(leftreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=*lpos+right[rightbase+j];
}
}
}
break;
case POW:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
const LSCALAR* lpos=left+(leftreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=pow(*lpos,right[rightbase+j]);
}
}
}
break;
case SUB:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
const LSCALAR* lpos=left+(leftreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=*lpos-right[rightbase+j];
}
}
}
break;
case MUL:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
const LSCALAR* lpos=left+(leftreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=*lpos*right[rightbase+j];
}
}
}
break;
case DIV:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
const LSCALAR* lpos=left+(leftreset?0:i*substep);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=*lpos/right[rightbase+j];
}
}
}
break;
default:
throw DataException("Unsupported binary operation");
}
}
template <>
void
binaryOpVectorLeftScalar(DataTypes::RealVectorType& res, // where result is to be stored
typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
const DataTypes::real_t* left, // LHS of calculation
const bool leftreset, // true if LHS is providing a single sample of 1 value only
const DataTypes::RealVectorType& right, // RHS of the calculation
typename DataTypes::RealVectorType::size_type rightOffset, // where to start reading RHS values
escript::ES_optype operation, // operation to perform
bool singlerightsample); // right consists of a single sample
template <class ResVEC, class LVEC, class RVEC>
void
binaryOpVector(ResVEC& res, // where result is to be stored
typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
const typename ResVEC::size_type sampleSize, // number of values in each sample
const LVEC& left, // LHS of calculation
typename LVEC::size_type leftOffset, // where to start reading LHS values
const bool leftreset, // Is LHS only supplying a single sample instead of a bunch of them
const RVEC& right, // RHS of the calculation
typename RVEC::size_type rightOffset, // where to start reading RHS values
const bool rightreset, // Is RHS only supplying a single sample instead of a bunch of them
escript::ES_optype operation) // operation to perform
{
switch (operation)
{
case ADD:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=left[leftbase+j]+right[rightbase+j];
}
}
}
break;
case POW:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],right[rightbase+j]);
}
}
}
break;
case SUB:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=left[leftbase+j]-right[rightbase+j];
}
}
}
break;
case MUL:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=left[leftbase+j]*right[rightbase+j];
}
}
}
break;
case DIV:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
{
typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
for (typename ResVEC::size_type j=0;j<sampleSize;++j)
{
res[i*sampleSize+resOffset+j]=left[leftbase+j]/right[rightbase+j];
}
}
}
break;
default:
throw DataException("Unsupported binary operation");
}
}
template <>
void
binaryOpVector(DataTypes::RealVectorType& res, // where result is to be stored
typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
const DataTypes::RealVectorType& left, // LHS of calculation
typename DataTypes::RealVectorType::size_type leftOffset, // where to start reading LHS values
const bool leftreset, // Is LHS only supplying a single sample instead of a bunch of them
const DataTypes::RealVectorType& right, // RHS of the calculation
typename DataTypes::RealVectorType::size_type rightOffset, // where to start reading RHS values
const bool rightreset, // Is RHS only supplying a single sample instead of a bunch of them
escript::ES_optype operation); // operation to perform
#define OPVECLAZYBODY(X) \
for (size_t j=0;j<onumsteps;++j)\
{\
for (size_t i=0;i<numsteps;++i,res+=resultStep) \
{ \
for (size_t s=0; s<chunksize; ++s)\
{\
res[s] = X;\
}\
/* tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X);*/ \
lroffset+=leftstep; \
rroffset+=rightstep; \
}\
lroffset+=oleftstep;\
rroffset+=orightstep;\
}
template <class ResELT, class LELT, class RELT>
void
binaryOpVectorLazyArithmeticHelper(ResELT* res,
const LELT* left,
const RELT* right,
const size_t chunksize,
const size_t onumsteps,
const size_t numsteps,
const size_t resultStep,
const size_t leftstep,
const size_t rightstep,
const size_t oleftstep,
const size_t orightstep,
size_t lroffset,
size_t rroffset,
escript::ES_optype operation) // operation to perform
{
switch (operation)
{
case ADD:
OPVECLAZYBODY((left[lroffset+s]+right[rroffset+s]));
break;
case POW:
OPVECLAZYBODY(pow(left[lroffset+s],right[rroffset+s]))
break;
case SUB:
OPVECLAZYBODY(left[lroffset+s]-right[rroffset+s])
break;
case MUL:
OPVECLAZYBODY(left[lroffset+s]*right[rroffset+s])
break;
case DIV:
OPVECLAZYBODY(left[lroffset+s]/right[rroffset+s])
break;
default:
ESYS_ASSERT(false, "Invalid operation. This should never happen!");
// I can't throw here because this will be called inside a parallel section
}
}
template <class ResELT, class LELT, class RELT>
void
binaryOpVectorLazyRelationalHelper(ResELT* res,
const LELT* left,
const RELT* right,
const size_t chunksize,
const size_t onumsteps,
const size_t numsteps,
const size_t resultStep,
const size_t leftstep,
const size_t rightstep,
const size_t oleftstep,
const size_t orightstep,
size_t lroffset,
size_t rroffset,
escript::ES_optype operation) // operation to perform
{
switch (operation)
{
case LESS:
OPVECLAZYBODY(left[lroffset+s]<right[rroffset+s])
break;
case GREATER:
OPVECLAZYBODY(left[lroffset+s]>right[rroffset+s])
break;
case GREATER_EQUAL:
OPVECLAZYBODY(left[lroffset+s]>=right[rroffset+s])
break;
case LESS_EQUAL:
OPVECLAZYBODY(left[lroffset+s]<=right[rroffset+s])
break;
default:
ESYS_ASSERT(false, "Invalid operation. This should never happen!");
// I can't throw here because this will be called inside a parallel section
}
}
/* trying to make a single version for all Tagged+Expanded interactions */
template <class ResVEC, class LVEC, class RVEC>
void
binaryOpVectorTagged(ResVEC& res, // where result is to be stored
const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
const typename ResVEC::size_type DPPSample, // number of datapoints per sample
const typename ResVEC::size_type DPSize, // datapoint size
const LVEC& left, // LHS of calculation
bool leftscalar,
const RVEC& right, // RHS of the calculation
bool rightscalar,
bool lefttagged, // true if left object is the tagged one
const DataTagged& tagsource, // where to get tag offsets from
escript::ES_optype operation) // operation to perform
{
typename ResVEC::size_type lstep=leftscalar?1:DPSize;
typename ResVEC::size_type rstep=rightscalar?1:DPSize;
typename ResVEC::size_type limit=samplesToProcess*DPPSample;
switch (operation)
{
case ADD:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<limit;++i)
{
typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
for (typename ResVEC::size_type j=0;j<DPSize;++j)
{
res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]+right[rightbase+j*(!rightscalar)];
}
}
}
break;
case POW:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<limit;++i)
{
typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
for (typename ResVEC::size_type j=0;j<DPSize;++j)
{
res[i*DPSize+j]=pow(left[leftbase+j*(!leftscalar)],right[rightbase+j*(!rightscalar)]);
}
}
}
break;
case SUB:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<limit;++i)
{
typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
for (typename ResVEC::size_type j=0;j<DPSize;++j)
{
res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]-right[rightbase+j*(!rightscalar)];
}
}
}
break;
case MUL:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<limit;++i)
{
typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
for (typename ResVEC::size_type j=0;j<DPSize;++j)
{
res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]*right[rightbase+j*(!rightscalar)];
}
}
}
break;
case DIV:
{
#pragma omp parallel for
for (typename ResVEC::size_type i=0;i<limit;++i)
{
typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
for (typename ResVEC::size_type j=0;j<DPSize;++j)
{
res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]/right[rightbase+j*(!rightscalar)];
}
}
}
break;
default:
throw DataException("Unsupported binary operation");
}
}
template<>
void
binaryOpVectorTagged(DataTypes::RealVectorType& res, // where result is to be stored
const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
const typename DataTypes::RealVectorType::size_type DPPSample, // number of datapoints per sample
const typename DataTypes::RealVectorType::size_type DPSize, // datapoint size
const DataTypes::RealVectorType& left, // LHS of calculation
const bool leftscalar,
const DataTypes::RealVectorType& right, // RHS of the calculation
const bool rightscalar,
const bool lefttagged, // true if left object is the tagged one
const DataTagged& tagsource, // where to get tag offsets from
escript::ES_optype operation);
template <class BinaryFunction>
inline
DataTypes::real_t
reductionOpVector(const DataTypes::RealVectorType& left,
const DataTypes::ShapeType& leftShape,
DataTypes::RealVectorType::size_type offset,
BinaryFunction operation,
DataTypes::real_t initial_value)
{
ESYS_ASSERT((left.size()>0)&&checkOffset(left,leftShape,offset),
"Couldn't perform reductionOp due to insufficient storage.");
DataTypes::real_t current_value=initial_value;
for (DataTypes::RealVectorType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
current_value=operation(current_value,left[offset+i]);
}
return current_value;
}
template <class BinaryFunction>
inline
DataTypes::real_t
reductionOpVector(const DataTypes::CplxVectorType& left,
const DataTypes::ShapeType& leftShape,
DataTypes::CplxVectorType::size_type offset,
BinaryFunction operation,
DataTypes::real_t initial_value)
{
ESYS_ASSERT((left.size()>0)&&checkOffset(left,leftShape,offset),
"Couldn't perform reductionOp due to insufficient storage.");
DataTypes::real_t current_value=initial_value;
for (DataTypes::RealVectorType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
current_value=operation(current_value,left[offset+i]);
}
return current_value;
}
int
matrix_inverse(const DataTypes::RealVectorType& in,
const DataTypes::ShapeType& inShape,
DataTypes::RealVectorType::size_type inOffset,
DataTypes::RealVectorType& out,
const DataTypes::ShapeType& outShape,
DataTypes::RealVectorType::size_type outOffset,
int count,
LapackInverseHelper& helper);
void
matrixInverseError(int err);
inline
bool
vectorHasNaN(const DataTypes::RealVectorType& in, DataTypes::RealVectorType::size_type inOffset, size_t count)
{
for (size_t z=inOffset;z<inOffset+count;++z)
{
if (nancheck(in[z]))
{
return true;
}
}
return false;
}
inline
bool
vectorHasNaN(const DataTypes::CplxVectorType& in, DataTypes::CplxVectorType::size_type inOffset, size_t count)
{
for (size_t z=inOffset;z<inOffset+count;++z)
{
if (nancheck(in[z]))
{
return true;
}
}
return false;
}
} // end namespace escript
#endif // __ESCRIPT_DATAMATHS_H__
escript::binaryOpVector
void binaryOpVector(DataTypes::RealVectorType &res, typename DataTypes::RealVectorType::size_type resOffset, const typename DataTypes::RealVectorType::size_type samplesToProcess, const typename DataTypes::RealVectorType::size_type sampleSize, const DataTypes::RealVectorType &left, typename DataTypes::RealVectorType::size_type leftOffset, const bool leftreset, const DataTypes::RealVectorType &right, typename DataTypes::RealVectorType::size_type rightOffset, const bool rightreset, escript::ES_optype operation)
Definition: DataVectorOps.cpp:767
escript::determineResultShape
DataTypes::ShapeType determineResultShape(const DataTypes::ShapeType &left, const DataTypes::ShapeType &right)
Determine the shape of the result array for a matrix multiplication of the given views.
Definition: DataVectorOps.cpp:169
escript::eigenvalues
void eigenvalues(const DataTypes::RealVectorType &in, const DataTypes::ShapeType &inShape, typename DataTypes::RealVectorType::size_type inOffset, DataTypes::RealVectorType &ev, const DataTypes::ShapeType &evShape, typename DataTypes::RealVectorType::size_type evOffset)
solves a local eigenvalue problem
Definition: DataVectorOps.h:638
escript::binaryOpVectorLazyRelationalHelper
void binaryOpVectorLazyRelationalHelper(ResELT *res, const LELT *left, const RELT *right, const size_t chunksize, const size_t onumsteps, const size_t numsteps, const size_t resultStep, const size_t leftstep, const size_t rightstep, const size_t oleftstep, const size_t orightstep, size_t lroffset, size_t rroffset, escript::ES_optype operation)
Definition: DataVectorOps.h:1219
escript::DataTypes::DataVectorAlt< real_t >::size_type
DataTypes::vec_size_type size_type
Definition: DataVectorAlt.h:49
escript::eigenvalues1
void eigenvalues1(const DataTypes::real_t A00, DataTypes::real_t *ev0)
solves a 1x1 eigenvalue A*V=ev*V problem
Definition: ArrayOps.h:160
escript::DataTagged::getPointOffset
virtual DataTypes::RealVectorType::size_type getPointOffset(int sampleNo, int dataPointNo) const
getPointOffset
Definition: DataTagged.cpp:897
escript
Definition: AbstractContinuousDomain.cpp:22
escript::transpose
void transpose(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset, int axis_offset)
Transpose each data point of this Data object around the given axis.
Definition: DataVectorOps.h:342
escript::ADD
Definition: ES_optype.h:30
escript::eigenvalues_and_eigenvectors3
void eigenvalues_and_eigenvectors3(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A02, const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *ev2, DataTypes::real_t *V00, DataTypes::real_t *V10, DataTypes::real_t *V20, DataTypes::real_t *V01, DataTypes::real_t *V11, DataTypes::real_t *V21, DataTypes::real_t *V02, DataTypes::real_t *V12, DataTypes::real_t *V22, const DataTypes::real_t tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition: ArrayOps.h:450
LapackInverseHelper.h
escript::vectorHasNaN
bool vectorHasNaN(const DataTypes::RealVectorType &in, DataTypes::RealVectorType::size_type inOffset, size_t count)
returns true if the vector contains NaN
Definition: DataVectorOps.h:1469
escript::LapackInverseHelper
Definition: LapackInverseHelper.h:26
escript::DataTypes::getRelIndex
vec_size_type getRelIndex(const DataTypes::ShapeType &shape, vec_size_type i)
Compute the offset (in 1D vector) of a given subscript with a shape.
Definition: DataTypes.h:233
escript::MUL
Definition: ES_optype.h:32
escript::binaryOpVectorLazyArithmeticHelper
void binaryOpVectorLazyArithmeticHelper(ResELT *res, const LELT *left, const RELT *right, const size_t chunksize, const size_t onumsteps, const size_t numsteps, const size_t resultStep, const size_t leftstep, const size_t rightstep, const size_t oleftstep, const size_t orightstep, size_t lroffset, size_t rroffset, escript::ES_optype operation)
Definition: DataVectorOps.h:1172
escript::DataTypes::ShapeType
std::vector< int > ShapeType
The shape of a single datapoint.
Definition: DataTypes.h:42
escript::binaryOpVectorLeftScalar
void binaryOpVectorLeftScalar(DataTypes::RealVectorType &res, typename DataTypes::RealVectorType::size_type resOffset, const typename DataTypes::RealVectorType::size_type samplesToProcess, const typename DataTypes::RealVectorType::size_type sampleSize, const DataTypes::real_t *left, const bool leftreset, const DataTypes::RealVectorType &right, typename DataTypes::RealVectorType::size_type rightOffset, escript::ES_optype operation, bool singlerightsample)
Definition: DataVectorOps.cpp:638
escript::DIV
Definition: ES_optype.h:33
escript::DataTagged
Simulates a full dataset accessible via sampleNo and dataPointNo.
Definition: DataTagged.h:44
DataTagged.h
escript::symmetric
void symmetric(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset)
computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
Definition: DataVectorOps.h:102
escript::DataTypes::DataVectorAlt::size
size_type size() const
Return the number of elements in this DataVectorAlt.
Definition: DataVectorAlt.h:214
escript::eigenvalues2
void eigenvalues2(const T A00, const T A01, const T A11, T *ev0, T *ev1)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
Definition: ArrayOps.h:182
DataException.h
escript::matMult
void matMult(const DataTypes::RealVectorType &left, const DataTypes::ShapeType &leftShape, DataTypes::RealVectorType::size_type leftOffset, const DataTypes::RealVectorType &right, const DataTypes::ShapeType &rightShape, DataTypes::RealVectorType::size_type rightOffset, DataTypes::RealVectorType &result, const DataTypes::ShapeType &resultShape)
Perform a matrix multiply of the given views.
Definition: DataVectorOps.cpp:38
ArrayOps.h
OPVECLAZYBODY
#define OPVECLAZYBODY(X)
Definition: DataVectorOps.h:1147
escript::hermitian
void hermitian(const DataTypes::CplxVectorType &in, const DataTypes::ShapeType &inShape, DataTypes::CplxVectorType::size_type inOffset, DataTypes::CplxVectorType &ev, const DataTypes::ShapeType &evShape, DataTypes::CplxVectorType::size_type evOffset)
computes an hermitian matrix from your square matrix A: (A + adjoint(A)) / 2
Definition: DataVectorOps.cpp:915
escript::binaryOpVectorRightScalar
void binaryOpVectorRightScalar(DataTypes::RealVectorType &res, typename DataTypes::RealVectorType::size_type resOffset, const typename DataTypes::RealVectorType::size_type samplesToProcess, const typename DataTypes::RealVectorType::size_type sampleSize, const DataTypes::RealVectorType &left, typename DataTypes::RealVectorType::size_type leftOffset, const DataTypes::real_t *right, const bool rightreset, escript::ES_optype operation, bool singleleftsample)
Definition: DataVectorOps.cpp:499
escript::DataTypes::DataVectorAlt< cplx_t >
escript::ES_optype
ES_optype
Definition: ES_optype.h:26
escript::matrix_inverse
int matrix_inverse(const DataTypes::RealVectorType &in, const DataTypes::ShapeType &inShape, DataTypes::RealVectorType::size_type inOffset, DataTypes::RealVectorType &out, const DataTypes::ShapeType &outShape, DataTypes::RealVectorType::size_type outOffset, int count, LapackInverseHelper &helper)
computes the inverses of square (up to 3x3) matricies
Definition: DataVectorOps.cpp:206
escript::binaryOpVectorTagged
void binaryOpVectorTagged(DataTypes::RealVectorType &res, const typename DataTypes::RealVectorType::size_type samplesToProcess, const typename DataTypes::RealVectorType::size_type DPPSample, const typename DataTypes::RealVectorType::size_type DPSize, const DataTypes::RealVectorType &left, const bool leftscalar, const DataTypes::RealVectorType &right, const bool rightscalar, const bool lefttagged, const DataTagged &tagsource, escript::ES_optype operation)
Definition: DataVectorOps.cpp:339
escript::eigenvalues_and_eigenvectors2
void eigenvalues_and_eigenvectors2(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A11, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *V00, DataTypes::real_t *V10, DataTypes::real_t *V01, DataTypes::real_t *V11, const DataTypes::real_t tol)
solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are ordered by increasing valu...
Definition: ArrayOps.h:341
escript::checkOffset
bool checkOffset(const VEC &data, const DataTypes::ShapeType &shape, typename VEC::size_type offset)
Definition: DataVectorOps.h:813
escript::eigenvalues_and_eigenvectors1
void eigenvalues_and_eigenvectors1(const DataTypes::real_t A00, DataTypes::real_t *ev0, DataTypes::real_t *V00, const DataTypes::real_t tol)
solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
Definition: ArrayOps.h:249
escript::GREATER_EQUAL
Definition: ES_optype.h:79
escript::LESS_EQUAL
Definition: ES_optype.h:80
escript::swapaxes
void swapaxes(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset, int axis0, int axis1)
swaps the components axis0 and axis1.
Definition: DataVectorOps.h:487
escript::matrixInverseError
void matrixInverseError(int err)
throws an appropriate exception based on failure of matrix_inverse.
Definition: DataVectorOps.cpp:185
V
#define V(_K_, _I_)
Definition: ShapeFunctions.cpp:120
escript::DataTypes::DataVectorAlt< real_t >::ElementType
real_t ElementType
Definition: DataVectorAlt.h:42
DataAbstract.h
escript::DataTypes::noValues
int noValues(const ShapeType &shape)
Calculate the number of values in a datapoint with the given shape.
Definition: DataTypes.cpp:90
escript::DataException
Definition: DataException.h:26
ESCRIPT_DLL_API
#define ESCRIPT_DLL_API
Definition: escriptcore/src/system_dep.h:29
escript::trace
void trace(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset, int axis_offset)
computes the trace of a matrix
Definition: DataVectorOps.h:241
escript::DataTypes::RealVectorType
escript::DataTypes::DataVectorAlt< real_t > RealVectorType
Vector to store underlying data.
Definition: DataVector.h:38
escript::reductionOpVector
DataTypes::real_t reductionOpVector(const DataTypes::RealVectorType &left, const DataTypes::ShapeType &leftShape, DataTypes::RealVectorType::size_type offset, BinaryFunction operation, DataTypes::real_t initial_value)
Perform the given data point reduction operation on the data point specified by the given offset into...
Definition: DataVectorOps.h:1393
escript::antisymmetric
void antisymmetric(const VEC &in, const DataTypes::ShapeType &inShape, typename VEC::size_type inOffset, VEC &ev, const DataTypes::ShapeType &evShape, typename VEC::size_type evOffset)
computes a antisymmetric matrix from your square matrix A: (A - transpose(A)) / 2 ...
Definition: DataVectorOps.h:151
escript::POW
Definition: ES_optype.h:34
escript::eigenvalues_and_eigenvectors
void eigenvalues_and_eigenvectors(const DataTypes::RealVectorType &in, const DataTypes::ShapeType &inShape, DataTypes::RealVectorType::size_type inOffset, DataTypes::RealVectorType &ev, const DataTypes::ShapeType &evShape, DataTypes::RealVectorType::size_type evOffset, DataTypes::RealVectorType &V, const DataTypes::ShapeType &VShape, DataTypes::RealVectorType::size_type VOffset, const double tol=1.e-13)
solves a local eigenvalue problem
Definition: DataVectorOps.h:746
escript::DataTypes::getRank
int getRank(const DataTypes::ShapeType &shape)
Return the rank (number of dimensions) of the given shape.
Definition: DataTypes.h:218
ESYS_ASSERT
#define ESYS_ASSERT(a, b)
EsysAssert is a MACRO that will throw an exception if the boolean condition specified is false...
Definition: Assert.h:78
escript::GREATER
Definition: ES_optype.h:78
escript::nancheck
bool nancheck(DataTypes::real_t d)
acts as a wrapper to isnan.
Definition: ArrayOps.h:116
escript::antihermitian
void antihermitian(const DataTypes::CplxVectorType &in, const DataTypes::ShapeType &inShape, typename DataTypes::CplxVectorType::size_type inOffset, DataTypes::CplxVectorType &ev, const DataTypes::ShapeType &evShape, typename DataTypes::CplxVectorType::size_type evOffset)
computes a antihermitian matrix from your square matrix A: (A - adjoint(A)) / 2
Definition: DataVectorOps.cpp:962
escript::SUB
Definition: ES_optype.h:31
escript::DataTypes::real_t
double real_t
type of all real-valued scalars in escript
Definition: DataTypes.h:50
escript::LESS
Definition: ES_optype.h:77
escript::eigenvalues3
void eigenvalues3(const DataTypes::real_t A00, const DataTypes::real_t A01, const DataTypes::real_t A02, const DataTypes::real_t A11, const DataTypes::real_t A12, const DataTypes::real_t A22, DataTypes::real_t *ev0, DataTypes::real_t *ev1, DataTypes::real_t *ev2)
solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
Definition: ArrayOps.h:206