1//
2// This file is part of Gambit
3// Copyright (c) 1994-2016, The Gambit Project (http://www.gambit-project.org)
4//
5// FILE: src/libgambit/dvector.imp
6// Implementation of doubly-partitioned vector class
7//
8// This program is free software; you can redistribute it and/or modify
9// it under the terms of the GNU General Public License as published by
10// the Free Software Foundation; either version 2 of the License, or
11// (at your option) any later version.
12//
13// This program is distributed in the hope that it will be useful,
14// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16// GNU General Public License for more details.
17//
18// You should have received a copy of the GNU General Public License
19// along with this program; if not, write to the Free Software
20// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
21//
22
23#include "dvector.h"
24
25namespace Gambit {
26
27//--------------------------------------------------------------------------
28//          DVector<T>: Private and protected member functions
29//--------------------------------------------------------------------------
30
31template <class T>
32int DVector<T>::sum(int part, const PVector<int> &v) const
33{
34  int s = 0;
35
36  Array<int> len(v.Lengths());
37  for (int j = 1; j <= len[part]; j++)
38    s += v(part, j);
39
40  return s;
41}
42
43template <class T> void DVector<T>::setindex(void)
44{
45  int index = 1;
46
47  for (int i = 1; i <= dvlen.Length(); i++)  {
48    dvptr[i] = this->svptr + index - 1;
49    dvidx[i] = index;
50    index += dvlen[i];
51  }
52}
53
54template <class T> bool DVector<T>::Check(const DVector<T> &v) const
55{
56  for (int i = 1; i <= dvlen.Length(); i++)
57    if (dvlen[i] != v.dvlen[i])   return false;
58  return true;
59}
60
61//--------------------------------------------------------------------------
62//    DVector<T>: Constructors, destructor, and constructive operators
63//--------------------------------------------------------------------------
64
65template <class T> DVector<T>::DVector(void) : dvptr(0)
66{ }
67
68template <class T> DVector<T>::DVector(const PVector<int> &sig)
69  : PVector<T>((Array<int>)sig), dvlen(sig.Lengths().Length()),
70    dvidx(sig.Lengths().Length())
71{
72  dvptr = new T **[dvlen.Length()];
73  dvptr -= 1;
74
75  for (int i = 1; i <= dvlen.Length(); i++)
76    dvlen[i] = sig.Lengths()[i];
77
78  setindex();
79}
80
81template <class T> DVector<T>::DVector(const Vector<T> &val,
82					   const PVector<int> &sig)
83  : PVector<T>(val, sig), dvlen(sig.Lengths().Length()),
84    dvidx(sig.Lengths().Length())
85{
86  dvptr = new T **[dvlen.Length()];
87  dvptr -= 1;
88
89  for (int i = 1; i <= dvlen.Length(); i++)
90    dvlen[i] = sig.Lengths()[i];
91
92  setindex();
93}
94
95template <class T> DVector<T>::DVector(const DVector<T> &v)
96  : PVector<T>(v), dvlen(v.dvlen), dvidx(v.dvidx)
97{
98  dvptr = new T **[dvlen.Length()];
99  dvptr -= 1;
100
101  setindex();
102}
103
104template <class T> DVector<T>::~DVector()
105{
106  if (dvptr)  delete [] (dvptr + 1);
107}
108
109template <class T> DVector<T> &DVector<T>::operator=(const DVector<T> &v)
110{
111  if (!Check(v)) {
112    throw DimensionException();
113  }
114
115  PVector<T>::operator=(v);
116  return *this;
117}
118
119template <class T> DVector<T> &DVector<T>::operator=(const PVector<T> &v)
120{
121  PVector<T>::operator=(v);
122  return *this;
123}
124
125template <class T> DVector<T> &DVector<T>::operator=(const Vector<T> &v)
126{
127  PVector<T>::operator=(v);
128  return *this;
129}
130
131template <class T> DVector<T> &DVector<T>::operator=(T c)
132{
133  PVector<T>::operator=(c);
134  return *this;
135}
136
137
138
139//--------------------------------------------------------------------------
140//                    DVector<T>: Operator definitions
141//--------------------------------------------------------------------------
142
143template <class T> T &DVector<T>::operator()(int a, int b, int c)
144{
145  if (dvlen.First() > a || a > dvlen.Last()) {
146    throw IndexException();
147  }
148  if (1 > b || b > dvlen[a]) {
149    throw IndexException();
150  }
151  if (1 > c || c > this->svlen[dvidx[a] + b - 1]) {
152    throw IndexException();
153  }
154  return dvptr[a][b][c];
155}
156
157template <class T> const T &DVector<T>::operator()(int a, int b, int c) const
158{
159  if (dvlen.First() > a || a > dvlen.Last()) {
160    throw IndexException();
161  }
162  if (1 > b || b > dvlen[a]) {
163    throw IndexException();
164  }
165  if (1 > c || c > this->svlen[dvidx[a] + b - 1]) {
166    throw IndexException();
167  }
168  return dvptr[a][b][c];
169}
170
171template <class T>
172DVector<T> DVector<T>::operator+(const DVector<T> &v) const
173{
174  if (!Check(v)) {
175    throw DimensionException();
176  }
177
178  DVector<T> tmp(*this);
179  tmp.PVector<T>::operator+=(v);
180  return tmp;
181}
182
183template <class T>
184DVector<T> &DVector<T>::operator+=(const DVector<T> &v)
185{
186  if (!Check(v)) {
187    throw DimensionException();
188  }
189  PVector<T>::operator+=(v);
190  return *this;
191}
192
193template <class T> DVector<T> DVector<T>::operator-(void) const
194{
195  DVector<T> tmp(*this);
196  for (int i = this->First(); i <= this->Last(); i++)
197    tmp[i] = -tmp[i];
198  return tmp;
199}
200
201template <class T>
202DVector<T> DVector<T>::operator-(const DVector<T> &v) const
203{
204  if (!Check(v)) {
205    throw DimensionException();
206  }
207  DVector<T> tmp(*this);
208  tmp.PVector<T>::operator-=(v);
209  return tmp;
210}
211
212template <class T>
213DVector<T> &DVector<T>::operator-=(const DVector<T> &v)
214{
215  if (!Check(v)) {
216    throw DimensionException();
217  }
218  PVector<T>::operator-=(v);
219  return *this;
220}
221
222template <class T> T DVector<T>::operator*(const DVector<T> &v) const
223{
224  if (!Check(v)) {
225    throw DimensionException();
226  }
227  return (*this).PVector<T>::operator*(v);
228}
229
230template <class T> DVector<T> &DVector<T>::operator*=(const T &c)
231{
232  PVector<T>::operator*=(c);
233  return *this;
234}
235
236template <class T> DVector<T> DVector<T>::operator/(const T &c) const
237{
238  DVector<T> tmp(*this);
239  tmp = tmp.PVector<T>::operator/(c);
240  return tmp;
241}
242
243template <class T> bool DVector<T>::operator==(const DVector<T> &v) const
244{
245  if (!Check(v)) {
246    throw DimensionException();
247  }
248  return PVector<T>::operator==(v);
249}
250
251template <class T> bool DVector<T>::operator!=(const DVector<T> &v) const
252{
253  return !(*this == v);
254}
255
256//-------------------------------------------------------------------------
257//                 DVector<T>: General data access
258//-------------------------------------------------------------------------
259
260template <class T>
261void DVector<T>::CopySubRow(int row, int col, const DVector<T> &v)
262{
263  if (!Check(v)) {
264    throw DimensionException();
265  }
266  if (dvlen.First() > row || row > dvlen.Last()) {
267    throw IndexException();
268  }
269  if (1 > col || col > dvlen[row]) {
270    throw IndexException();
271  }
272
273  for (int i = 1; i <= this->svlen[dvidx[row]+col-1]; i++)
274    dvptr[row][col][i] = v.dvptr[row][col][i];
275}
276
277
278
279template <class T> const Array<int> &DVector<T>::DPLengths(void) const
280{
281  return dvlen;
282}
283
284} // end namespace Gambit
285