1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
/***************************************************************************
* blitz/matutri.h Declarations for UpperTriangular matrices
*
* $Id: matutri.h 1414 2005-11-01 22:04:59Z cookedm $
*
* Copyright (C) 1997-2001 Todd Veldhuizen <tveldhui@oonumerics.org>
*
* This code was relicensed under the modified BSD license for use in SciPy
* by Todd Veldhuizen (see LICENSE.txt in the weave directory).
*
*
* Suggestions: blitz-dev@oonumerics.org
* Bugs: blitz-bugs@oonumerics.org
*
* For more information, please see the Blitz++ Home Page:
* http://oonumerics.org/blitz/
*
***************************************************************************/
#ifndef BZ_MATUTRI_H
#define BZ_MATUTRI_H
#ifndef BZ_MSTRUCT_H
#error <blitz/matutri.h> must be included via <blitz/mstruct.h>
#endif
BZ_NAMESPACE(blitz)
// Upper triangular, column major ordering
// [ 0 1 3 6 ]
// [ . 2 4 7 ]
// [ . . 5 8 ]
// [ . . . 9 ]
class UpperTriangularIterator {
public:
UpperTriangularIterator(unsigned rows, unsigned cols)
{
BZPRECONDITION(rows == cols);
size_ = rows;
good_ = true;
offset_ = 0;
i_ = 0;
j_ = 0;
}
operator bool() const { return good_; }
void operator++()
{
BZPRECONDITION(good_);
++offset_;
++i_;
if (i_ > j_)
{
i_ = 0;
++j_;
if (j_ == size_)
good_ = false;
}
}
unsigned row() const
{ return i_; }
unsigned col() const
{ return j_; }
unsigned offset() const
{ return offset_; }
protected:
unsigned size_;
bool good_;
unsigned offset_;
unsigned i_, j_;
};
class UpperTriangular : public MatrixStructure {
public:
typedef UpperTriangularIterator T_iterator;
UpperTriangular()
: size_(0)
{ }
UpperTriangular(unsigned size)
: size_(size)
{ }
UpperTriangular(unsigned rows, unsigned cols)
: size_(rows)
{
BZPRECONDITION(rows == cols);
}
unsigned columns() const
{ return size_; }
unsigned coordToOffset(unsigned i, unsigned j) const
{
BZPRECONDITION(inRange(i,j));
BZPRECONDITION(j >= i);
return j*(j+1)/2 + i;
}
unsigned firstInRow(unsigned i) const
{ return 0; }
template<typename T_numtype>
T_numtype get(const T_numtype * restrict data,
unsigned i, unsigned j) const
{
BZPRECONDITION(inRange(i,j));
if (j >= i)
return data[coordToOffset(i,j)];
else
return ZeroElement<T_numtype>::zero();
}
template<typename T_numtype>
T_numtype& get(T_numtype * restrict data, unsigned i, unsigned j)
{
BZPRECONDITION(inRange(i,j));
if (j >= i)
return data[coordToOffset(i,j)];
else
return ZeroElement<T_numtype>::zero();
}
unsigned lastInRow(unsigned i) const
{ return size_ - 1; }
unsigned firstInCol(unsigned j) const
{ return 0; }
unsigned lastInCol(unsigned j) const
{ return j; }
bool inRange(const unsigned i,const unsigned j) const { return (i<size_) && (j<size_); }
unsigned numElements() const { return size_ * (size_ + 1) / 2; }
unsigned rows() const { return size_; }
void resize(const unsigned size) { size_ = size; }
void resize(const unsigned rows,const unsigned cols) {
BZPRECONDITION(rows == cols);
size_ = rows;
}
private:
unsigned size_;
};
BZ_NAMESPACE_END
#endif // BZ_MATUTRI_H
|
