← Back to branch summary

~durga/maus/online

« back to all changes in this revision

Viewing changes to tests/integration/optics/src/Tensor.cc

Builds but segv on executation of test_cpp

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/integration/optics/src/Tensor.cc
1
 
#include "Tensor.hh"
2
 
#include "Tensor3.hh"
3
 
 
4
 
const Squeal Tensor::_outOfRange = Squeal(Squeal::recoverable, "Tensor index out of range",       "Tensor::Tensor");
5
 
const Squeal Tensor::_outOfRank  = Squeal(Squeal::recoverable, "Tensor rank out of range",        "Tensor::Tensor");
6
 
const Squeal Tensor::_lowRank    = Squeal(Squeal::recoverable, "Do not use Tensor for rank <= 3", "Tensor::Tensor");
7
 
 
8
 
 
9
 
Tensor::Tensor(std::vector<int> size, double val) : _data(NULL), _size(size)
10
 
{
11
 
        SetTensor(size, val);
12
 
}
13
 
 
14
 
Tensor::Tensor(int size[],  int rank, double val) : _data(NULL)
15
 
{
16
 
        std::vector<int> sizeV(size, size+rank);
17
 
        SetTensor(sizeV, val);
18
 
}
19
 
 
20
 
Tensor::Tensor(int size,    int rank, double val) : _data(NULL)
21
 
{
22
 
        std::vector<int> sizeV(rank, size);
23
 
        SetTensor(sizeV, val);
24
 
}
25
 
 
26
 
Tensor::~Tensor()
27
 
{
28
 
        for(int i=0; i<_size[0]; i++)
29
 
                if(_data[i]) delete _data[i];
30
 
        if(_data && _size[0]>0) delete [] _data;
31
 
}
32
 
 
33
 
void Tensor::SetTensor(std::vector<int> size, double val)
34
 
{
35
 
        _size = size;
36
 
        if(size.size() < 4) throw(_lowRank);
37
 
        else if(size.size() == 4)
38
 
        {
39
 
                _data = new Tensor*[size[0]];
40
 
                for(int i=0; i<size[0]; i++)
41
 
                        _data[i] = new Tensor3(size[1], size[2], size[3], val);
42
 
        }
43
 
        else
44
 
        {
45
 
                _data = new Tensor*[size[0]];
46
 
                std::vector<int> newSize;
47
 
                for(unsigned int i=1; i<size.size(); i++)
48
 
                        newSize.push_back(size[i]);
49
 
                for(int i=0; i<size[0]; i++)
50
 
                        _data[i] = new Tensor(newSize, val);
51
 
        }
52
 
 
53
 
        SetAllFields(val);
54
 
        _size   = size;
55
 
}
56
 
 
57
 
void Tensor::SetAllFields(double value)
58
 
{
59
 
        for(int i=0; i<_size[0]; i++)
60
 
                _data[i]->SetAllFields(value);
61
 
 
62
 
}
63
 
 
64
 
void Tensor::Print(std::ostream& out) const
65
 
{
66
 
        for(unsigned int i=0; i<_size.size(); i++)
67
 
                out << "\n";
68
 
        for(int i=0; i<_size[0]; i++)
69
 
                out << (*_data[i]);
70
 
}
71
 
 
72
 
std::ostream& operator<<(std::ostream& out, const Tensor& t1)
73
 
{
74
 
        t1.Print(out); 
75
 
        return out;
76
 
}
77
 
 
78
 
bool Tensor::IsZero() const
79
 
{
80
 
        for(int i=0; i<_size[0]; i++)
81
 
                if(!_data[i]->IsZero()) return false;
82
 
        return true;
83
 
}
84
 
 
85
 
Tensor operator/(const Tensor& t1,    const double& d1)
86
 
{
87
 
        Tensor tOut = t1;
88
 
        for(int i=0; i<t1._size[0]; i++) 
89
 
                (*tOut._data[i]) = (*tOut._data[i])/d1; 
90
 
        return tOut;
91
 
}
92
 
 
93
 
Tensor Tensor::Poisson(Tensor t1) const
94
 
{
95
 
        //:T_i::T_j: is of order i+j-3
96
 
        int    rankOut = GetRank() + t1.GetRank() - 3;
97
 
        Tensor tOut(6, rankOut, 0);
98
 
        //:T_i::T_j: = dT_i/du_p d/dConj{u_p} [dT_j/du_q d/dConj{u_q}]
99
 
        //so tOut[i1][i2]...[ik] = dGdUi(j, v)*Delta(i, GetConjugate(j))*pow(-1, j+1);
100
 
        
101
 
        //THIS IS THE BIT I AM WORKING ON - TRYING TO FIGURE IT OUT!
102
 
 
103
 
        return tOut;
104
 
}
105
 
 
106
 
/*
107
 
std::vector<Tensor*> Tensor::Exponential(int maxOrder) const
108
 
{
109
 
        //:T_i:^n is of order n*(i-2)+1
110
 
        //maxPower can be non-integer
111
 
        double  maxPower  = (maxOrder-1)/(double)(GetRank()-2);
112
 
        int     nFact     = 1; //n!
113
 
        std::vector<Tensor*> out; // sum(:t:^n)
114
 
        out.push_back(this->Clone());
115
 
        for(int n=2; (double)n<=maxPower; n++)
116
 
        {
117
 
                nFact *= n;
118
 
                Tensor* current = out.back(); 
119
 
                out.push_back(new Tensor(current->Poisson(*this)/nFact));
120
 
        }       
121
 
        return out;
122
 
}
123
 
 
124
 
CLHEP::HepVector Tensor::ExponentialPoisson(CLHEP::HepVector v, int maxOrder) const
125
 
{
126
 
        std::vector<Tensor*> ExpT = Exponential(maxOrder);
127
 
        CLHEP::HepVector out;
128
 
        for(unsigned int i=0; i<ExpT.size(); i++)
129
 
                out += ExpT[i]->Poisson(v);
130
 
        return out;
131
 
}
132
 
 
133
 
CLHEP::HepVector Tensor::Poisson(CLHEP::HepVector v) const;
134
 
{
135
 
        int              dimension = v.num_row();
136
 
        CLHEP::HepVector vOut(dimension);
137
 
        std::vector<int> index(_size.size());
138
 
        return vOut;
139
 
}
140
 
*/
141
 
/*
142
 
        CLHEP::HepVector     Poisson(CLHEP::HepVector v) const;
143
 
        CLHEP::HepMatrix     Poisson(CLHEP::HepMatrix m) const;
144
 
        Tensor               Poisson(Tensor           t) const;
145
 
*/
146
 
 
147
 
//Turn int index into an index for element of symmetric tensor
148
 
std::vector<int> Tensor::IndexToIndices(int index, int tensorSize, int startDimension)
149
 
{
150
 
        index++;
151
 
        std::vector<int> indices(1,-1);
152
 
        std::vector<int> subIndices(1,0);
153
 
        int hex        = tensorSize ;
154
 
        for(int i=0; i<index; i++)
155
 
        {
156
 
//              for(unsigned int i=0; i<indices.size(); i++) std::cout << indices[i] << " ";
157
 
                if     (indices.front() == hex-1) { indices = std::vector<int>(indices.size()+1, 0); indices.back()--;}
158
 
                if     (indices.back()  == hex-1)
159
 
                {
160
 
                        int j=indices.size()-1;
161
 
                        while(indices[j] == hex-1) j--;
162
 
                        for(int k=indices.size()-1; k>=j; k--) indices[k] = indices[j]+1;
163
 
                }
164
 
                else    indices.back()++;
165
 
        }
166
 
//      std::cout << std::endl;
167
 
        for(unsigned int i=0; i<indices.size(); i++) indices[i] += startDimension;
168
 
        return indices;
169
 
}
170
 
/*
171
 
        std::vector<int> indices(0);
172
 
        //convert from decimal to heximal
173
 
        int hex        = tensorSize - 1;
174
 
        for(int i=0; i<=index; i++)
175
 
        {
176
 
                int j = indices.size()-1;
177
 
                if(j >= 0) while(indices[j]  == hex && j>=0) j--;
178
 
                if(j >= 0) indices[j]++;
179
 
                else       indices = std::vector<int>(indices.size()+1);
180
 
        }
181
 
        return indices;
182
 
}
183
 
*/
184
 
 
185
 
double Tensor::Get(std::vector<int> position) const
186
 
{
187
 
        if(int(position.size()) != GetRank()) throw(Squeal(Squeal::recoverable, "Get rank different to Tensor rank", "Tensor::Get"));
188
 
        const Tensor* val = this;
189
 
        if(GetRank()>3)
190
 
        {
191
 
                for(unsigned int i=position.size()-1; i>3; i--)
192
 
                        val = val->_data[position[i]];
193
 
                val = val->_data[position[3]];
194
 
        }
195
 
        return ( *((Tensor3*)val) )[position[0]][position[1]][position[2]];
196
 
}
197
 
 
198
 
void   Tensor::Set     (std::vector<int> position, double value)
199
 
{
200
 
        Tensor* val = this;
201
 
        if(GetRank()>3)
202
 
        {
203
 
                for(unsigned int i=position.size()-1; i>3; i--)
204
 
                        val = val->_data[position[i]];
205
 
                val = val->_data[position[3]];
206
 
        }
207
 
        ( *((Tensor3*)val) )[position[0]][position[1]][position[2]] = value;
208
 
}
209
 
 
210