Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F67143875
SparseMatrix-inl.h
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Jun 20, 06:01
Size
33 KB
Mime Type
text/x-c++
Expires
Sat, Jun 22, 06:01 (2 d)
Engine
blob
Format
Raw Data
Handle
18309986
Attached To
rLAMMPS lammps
SparseMatrix-inl.h
View Options
#ifndef SPARSEMATRIX_INL_H
#define SPARSEMATRIX_INL_H
template
<
typename
T
>
TRI_COORD
<
T
>::
TRI_COORD
(
unsigned
row
,
unsigned
col
)
:
i
(
row
),
j
(
col
)
{}
template
<
typename
T
>
TRI_COORD
<
T
>::
TRI_COORD
(
unsigned
row
,
unsigned
col
,
T
val
,
bool
add_to
)
:
i
(
row
),
j
(
col
),
v
(
val
),
add
(
add_to
)
{}
// General flat index by value operator (by nth nonzero)
template
<
typename
T
>
inline
T
SparseMatrix
<
T
>::
operator
[](
INDEX
i
)
const
{
VICK
(
i
);
return
_val
[
i
];
}
// General flat index by reference operator (by nth nonzero)
template
<
typename
T
>
inline
T
&
SparseMatrix
<
T
>::
operator
[](
INDEX
i
)
{
VICK
(
i
);
return
_val
[
i
];
}
template
<
typename
T
>
T
SparseMatrix
<
T
>::
_zero
=
T
(
0
);
//-----------------------------------------------------------------------------
// triplet comparison operator returns true if x < y
//-----------------------------------------------------------------------------
template
<
typename
T
>
bool
triplet_comparision
(
const
TRI_COORD
<
T
>
&
x
,
const
TRI_COORD
<
T
>
&
y
)
{
const
bool
row_less
=
(
x
.
i
)
<
(
y
.
i
);
const
bool
row_equal
=
(
x
.
i
)
==
(
y
.
i
);
const
bool
col_less
=
(
x
.
j
)
<
(
y
.
j
);
return
(
row_less
||
(
row_equal
&&
col_less
));
}
//-----------------------------------------------------------------------------
// triplet comparison operator returns true if x == y
//-----------------------------------------------------------------------------
template
<
typename
T
>
bool
triplets_equal
(
const
TRI_COORD
<
T
>
&
x
,
const
TRI_COORD
<
T
>
&
y
)
{
return
x
.
i
==
y
.
i
&&
x
.
j
==
y
.
j
;
}
//-----------------------------------------------------------------------------
// multiply sparse matrix by a vector
//-----------------------------------------------------------------------------
template
<
typename
T
>
DenseVector
<
T
>
operator
*
(
const
SparseMatrix
<
T
>
&
A
,
const
Vector
<
T
>&
x
)
{
SparseMatrix
<
T
>::
compress
(
A
);
GCK
(
A
,
x
,
A
.
nCols
()
!=
x
.
size
(),
"SparseMatrix * Vector"
)
DenseVector
<
T
>
y
(
A
.
nRows
(),
true
);
INDEX
i
,
j
;
for
(
i
=
0
;
i
<
A
.
_nRowsCRS
;
i
++
)
for
(
j
=
A
.
_ia
[
i
];
j
<
A
.
_ia
[
i
+
1
];
j
++
)
y
(
i
)
+=
A
.
_val
[
j
]
*
x
(
A
.
_ja
[
j
]);
return
y
;
}
//-----------------------------------------------------------------------------
// multiply a vector by a sparse matrix
//-----------------------------------------------------------------------------
template
<
typename
T
>
DenseVector
<
T
>
operator
*
(
const
Vector
<
T
>&
x
,
const
SparseMatrix
<
T
>
&
A
)
{
return
A
.
transMat
(
x
);
}
//-----------------------------------------------------------------------------
// multiply sparse matrix by dense matrix
//-----------------------------------------------------------------------------
template
<
typename
T
>
DenseMatrix
<
T
>
operator
*
(
const
SparseMatrix
<
T
>
&
A
,
const
Matrix
<
T
>&
D
)
{
GCK
(
A
,
D
,
A
.
nCols
()
!=
D
.
nRows
(),
"SparseMatrix * DenseMatrix"
)
DenseMatrix
<
T
>
C
(
A
.
nRows
(),
D
.
nCols
(),
true
);
// initialized to zero
const
INDEX
J
=
D
.
nCols
();
INDEX
i
,
ik
,
j
,
k
;
for
(
i
=
0
;
i
<
A
.
_nRowsCRS
;
i
++
)
for
(
ik
=
A
.
_ia
[
i
];
ik
<
A
.
_ia
[
i
+
1
];
ik
++
)
for
(
j
=
0
;
j
<
J
;
j
++
)
C
(
i
,
j
)
+=
A
.
_val
[
ik
]
*
D
(
A
.
_ja
[
ik
],
j
);
// C(i,j) = S(i,k) * D(k, j)
return
C
;
}
//-----------------------------------------------------------------------------
// multiply sparse matrix by a diagonal matrix - scales each column
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>
operator
*
(
const
SparseMatrix
<
T
>
&
A
,
const
DiagonalMatrix
<
T
>&
D
)
{
GCK
(
A
,
D
,
A
.
nCols
()
!=
D
.
nRows
(),
"SparseMatrix * DiagonalMatrix"
)
SparseMatrix
<
T
>
C
(
A
);
// C has same sparcity as A
// C(i,j) = A(i,k) * D(k, j) * j==k
INDEX
i
,
ij
,
j
;
for
(
i
=
0
;
i
<
A
.
_nRowsCRS
;
i
++
)
for
(
ij
=
A
.
_ia
[
i
];
ij
<
A
.
_ia
[
i
+
1
];
ij
++
)
C
[
ij
]
=
A
.
_val
[
ij
]
*
D
(
A
.
_ja
[
ij
],
A
.
_ja
[
ij
]);
return
C
;
}
//-----------------------------------------------------------------------------
// multiplies two sparse matrices - assumes their output is sparse
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>
operator
*
(
const
SparseMatrix
<
T
>
&
A
,
const
SparseMatrix
<
T
>
&
B
)
{
SparseMatrix
<
T
>
At
(
A
.
transpose
());
SparseMatrix
<
T
>::
compress
(
B
);
GCK
(
A
,
B
,
A
.
nCols
()
!=
B
.
nRows
(),
"SparseMatrix * SparseMatrix"
);
SparseMatrix
<
T
>
C
(
A
.
nRows
(),
B
.
nCols
());
if
(
At
.
empty
()
||
B
.
empty
())
return
C
;
INDEX
k
,
ki
,
kj
;
INDEX
K
=
std
::
min
(
At
.
_nRowsCRS
,
B
.
_nRowsCRS
);
for
(
k
=
0
;
k
<
K
;
k
++
)
// loop over rows of A or B (smallest)
for
(
ki
=
At
.
_ia
[
k
];
ki
<
At
.
_ia
[
k
+
1
];
ki
++
)
// loop over row nonzeros of A
for
(
kj
=
B
.
_ia
[
k
];
kj
<
B
.
_ia
[
k
+
1
];
kj
++
)
// loop over row nonzeros of B
C
.
add
(
At
.
_ja
[
ki
],
B
.
_ja
[
kj
],
At
[
ki
]
*
B
[
kj
]);
// C(i,j) = At(k,i)*B(k, j)
C
.
compress
();
return
C
;
}
//-----------------------------------------------------------------------------
// default constructor - creates an empty sparsematrix with specified size
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>::
SparseMatrix
(
INDEX
rows
,
INDEX
cols
)
:
_val
(
NULL
),
_ia
(
NULL
),
_ja
(
NULL
),
_size
(
0
),
_nRowsCRS
(
0
),
_nRows
(
rows
),
_nCols
(
cols
)
{}
//-----------------------------------------------------------------------------
// copy constructor
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>::
SparseMatrix
(
const
SparseMatrix
<
T
>&
C
)
:
_val
(
NULL
),
_ia
(
NULL
),
_ja
(
NULL
)
{
_copy
(
C
);
}
//-----------------------------------------------------------------------------
// copy constructor - converts from DenseMatrix
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>::
SparseMatrix
(
const
DenseMatrix
<
T
>&
C
)
:
_val
(
NULL
),
_ia
(
NULL
),
_ja
(
NULL
)
{
reset
(
C
);
}
//-----------------------------------------------------------------------------
// destructor, cleans up internal storage
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>::~
SparseMatrix
()
{
_delete
();
}
//-----------------------------------------------------------------------------
// assigns internal storage for CRS
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
_create
(
INDEX
size
,
INDEX
nrows
)
{
_size
=
size
;
_nRowsCRS
=
nrows
;
// assign memory to hold matrix
try
{
_val
=
(
_size
*
nrows
)
?
new
T
[
_size
]
:
NULL
;
_ia
=
(
_size
*
nrows
)
?
new
INDEX
[
_nRowsCRS
+
1
]
:
NULL
;
_ja
=
(
_size
*
nrows
)
?
new
INDEX
[
_size
]
:
NULL
;
}
catch
(
std
::
exception
&
e
)
{
cout
<<
"Could not allocate SparseMatrix of "
<<
_size
<<
" nonzeros.
\n
"
;
exit
(
EXIT_FAILURE
);
}
if
(
!
_ia
)
return
;
// automatically handle the ends of rowpointer
*
_ia
=
0
;
// first non-zero is the zero index
_ia
[
_nRowsCRS
]
=
_size
;
// last row pointer is the size
}
//-----------------------------------------------------------------------------
// cleans up internal storage, but retains nRows & nCols
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
_delete
()
{
vector
<
TRI_COORD
<
T
>
>
().
swap
(
_tri
);
// completely deletes _tri
delete
[]
_val
;
delete
[]
_ia
;
delete
[]
_ja
;
_size
=
_nRowsCRS
=
0
;
_val
=
NULL
;
_ia
=
_ja
=
NULL
;
}
//-----------------------------------------------------------------------------
// full memory copy of C into this
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
_copy
(
const
SparseMatrix
<
T
>
&
C
)
{
compress
(
C
);
_delete
();
_create
(
C
.
size
(),
C
.
_nRowsCRS
);
if
(
_size
)
{
std
::
copy
(
C
.
_val
,
C
.
_val
+
_size
,
_val
);
std
::
copy
(
C
.
_ja
,
C
.
_ja
+
_size
,
_ja
);
}
if
(
_nRowsCRS
)
{
std
::
copy
(
C
.
_ia
,
C
.
_ia
+
_nRowsCRS
+
1
,
_ia
);
}
_nCols
=
C
.
_nCols
;
_nRows
=
C
.
_nRows
;
}
//----------------------------------------------------------------------------
// general sparse matrix assignment
//----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
_set_equal
(
const
Matrix
<
T
>
&
r
)
{
this
->
resize
(
r
.
nRows
(),
r
.
nCols
());
const
Matrix
<
T
>
*
ptr_r
=
&
r
;
const
SparseMatrix
<
T
>
*
s_ptr
=
dynamic_cast
<
const
SparseMatrix
<
T
>*>
(
ptr_r
);
if
(
s_ptr
)
this
->
reset
(
*
s_ptr
);
else
if
(
dynamic_cast
<
const
DiagonalMatrix
<
T
>*>
(
ptr_r
))
for
(
INDEX
i
=
0
;
i
<
r
.
size
();
i
++
)
set
(
i
,
i
,
r
[
i
]);
else
if
(
dynamic_cast
<
const
DenseMatrix
<
T
>*>
(
ptr_r
))
this
->
reset
(
r
);
else
{
cout
<<
"Error in general sparse matrix assignment
\n
"
;
exit
(
1
);
}
}
//-----------------------------------------------------------------------------
// returns the first row number with a nonzero entry or -1 if no rows
//-----------------------------------------------------------------------------
template
<
typename
T
>
int
SparseMatrix
<
T
>::
_first_nonzero_row_crs
()
const
{
if
(
!
_nRowsCRS
)
return
-
1
;
INDEX
r
;
for
(
r
=
0
;
r
<
_nRowsCRS
;
r
++
)
if
(
_ia
[
r
+
1
]
>
0
)
return
r
;
return
-
1
;
}
//-----------------------------------------------------------------------------
// converts T to CRS
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
compress
(
const
SparseMatrix
<
T
>
&
C
)
{
const_cast
<
SparseMatrix
<
T
>*>
(
&
C
)
->
compress
();
}
//-----------------------------------------------------------------------------
// merges all the _tri triples with CRS storage
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
compress
()
{
if
(
_tri
.
empty
())
return
;
// sort and find the number of unique triplets
const
INDEX
nUnique
=
CountUniqueTriplets
();
// max number of rows in new CRS structure
const
INDEX
nRows
=
std
::
max
((
INDEX
)
_tri
.
back
().
i
+
1
,
_nRowsCRS
);
// make a new CRS structure
INDEX
*
ia
=
new
INDEX
[
nRows
+
1
];
INDEX
*
ja
=
new
INDEX
[
nUnique
];
T
*
val
=
new
T
[
nUnique
];
ia
[
0
]
=
0
;
INDEX
i
;
for
(
i
=
1
;
i
<
nRows
;
i
++
)
ia
[
i
]
=~
0
;
// ~0 is max(INDEX)
ia
[
nRows
]
=
nUnique
;
INDEX
crs_pt
,
crs_row
,
tri_ct
;
// get the first CRS and triplet coordinates (if they exist)
TRI_COORD
<
T
>
nextCRS
,
nextTRI
(
_tri
[
0
]),
next
;
int
first_row
=
_first_nonzero_row_crs
();
if
(
first_row
!=
-
1
)
nextCRS
=
TRI_COORD
<
T
>
(
first_row
,
_ja
[
0
],
_val
[
0
]);
// merge sorted triplets into a new CRS structure
crs_pt
=
crs_row
=
tri_ct
=
0
;
// initialize counters
for
(
i
=
0
;
i
<
nUnique
;
i
++
)
{
// is the next non-zero in the new triplet vector
if
(
tri_ct
<
_tri
.
size
()
&&
(
triplet_comparision
(
nextTRI
,
nextCRS
)
||
crs_pt
>=
_size
))
{
next
=
nextTRI
;
// advance the triplet counter, and skip voided TRIPLET entries
do
tri_ct
++
;
while
(
tri_ct
<
_tri
.
size
()
&&
!~
_tri
[
tri_ct
].
j
);
// if not at the end of the vector, set the next triplet
if
(
tri_ct
<
_tri
.
size
())
nextTRI
=
_tri
[
tri_ct
];
}
// is the next nonzero in the old CRS data
else
if
(
crs_pt
<
_size
)
{
next
=
nextCRS
;
// advance the CRS counter, skip if it was the last one
if
(
++
crs_pt
>=
_size
)
continue
;
crs_row
+=
_ia
[
crs_row
+
1
]
==
crs_pt
;
// did the row advance
nextCRS
=
TRI_COORD
<
T
>
(
crs_row
,
_ja
[
crs_pt
],
_val
[
crs_pt
]);
}
else
cout
<<
"SparseMatrix - Error in compressing CRS
\n
"
;
// add next to the new CRS structure
// is this a new row (is j>0 and is ja[j] == 0)?
if
(
ia
[
next
.
i
]
==~
0
)
ia
[
next
.
i
]
=
i
;
ja
[
i
]
=
next
.
j
;
val
[
i
]
=
next
.
v
;
}
// sweep backwards through row pointers and check for skipped rows
for
(
i
=
nRows
-
1
;
i
>
0
;
i
--
)
ia
[
i
]
=
(
ia
[
i
]
==~
0
)
?
ia
[
i
+
1
]
:
ia
[
i
];
_delete
();
_val
=
val
;
_ia
=
ia
;
_ja
=
ja
;
_size
=
nUnique
;
_nRowsCRS
=
nRows
;
}
//-----------------------------------------------------------------------------
// Sorts the triplets, condenses duplicates, and returns the # of unique values
//-----------------------------------------------------------------------------
template
<
typename
T
>
INDEX
SparseMatrix
<
T
>::
CountUniqueTriplets
()
{
std
::
sort
(
_tri
.
begin
(),
_tri
.
end
(),
triplet_comparision
<
T
>
);
INDEX
i
,
nUnique
=
1
+
_size
;
// count unique entries
for
(
i
=
_tri
.
size
()
-
1
;
i
>
0
;
i
--
)
{
// for all new triplets
if
(
triplets_equal
(
_tri
[
i
-
1
],
_tri
[
i
]))
{
// previous has same index?
if
(
_tri
[
i
].
add
)
_tri
[
i
-
1
].
v
+=
_tri
[
i
].
v
;
// add to previous
else
_tri
[
i
-
1
].
v
=
_tri
[
i
].
v
;
// replace previous
_tri
[
i
].
j
=
~
0
;
// void this entry's col
}
else
nUnique
++
;
}
return
nUnique
;
}
//-----------------------------------------------------------------------------
// Index by copy operator - return zero if not found
//-----------------------------------------------------------------------------
template
<
typename
T
>
T
SparseMatrix
<
T
>::
operator
()(
INDEX
i
,
INDEX
j
)
const
{
MICK
(
i
,
j
);
// Matrix Index ChecKing
compress
(
*
this
);
if
(
i
>=
_nRowsCRS
||
_ia
[
i
+
1
]
==
_ia
[
i
])
return
0.0
;
unsigned
f
=
std
::
lower_bound
(
_ja
+
_ia
[
i
],
_ja
+
_ia
[
i
+
1
]
-
1
,
j
)
-
_ja
;
if
(
f
>=
_ia
[
i
]
&&
f
<
_ia
[
i
+
1
]
&&
_ja
[
f
]
==
j
)
return
_val
[
f
];
return
0.0
;
}
//-----------------------------------------------------------------------------
// Index by reference operator - add to _tri if not found
//-----------------------------------------------------------------------------
template
<
typename
T
>
T
&
SparseMatrix
<
T
>::
operator
()(
INDEX
i
,
INDEX
j
)
{
compress
(
*
this
);
MICK
(
i
,
j
);
// Matrix Index ChecKing
if
(
i
<
_nRowsCRS
&&
_ia
[
i
+
1
]
>
_ia
[
i
])
{
unsigned
f
=
std
::
lower_bound
(
_ja
+
_ia
[
i
],
_ja
+
_ia
[
i
+
1
]
-
1
,
j
)
-
_ja
;
if
(
f
>=
_ia
[
i
]
&&
f
<
_ia
[
i
+
1
]
&&
_ja
[
f
]
==
j
)
return
_val
[
f
];
}
// NEVER use index operator as LHS to modify values not already in the
// sparcity pattern - the crude check below will only catch this on the
// second infraction.
if
(
_zero
!=
T
(
0
))
cout
<<
"Use add or set for SparseMatrix
\n
"
;
return
_zero
;
}
//-----------------------------------------------------------------------------
// Sets (i,j) to value
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
set
(
INDEX
i
,
INDEX
j
,
T
v
)
{
MICK
(
i
,
j
);
// Matrix Index ChecKing
if
(
i
<
_nRowsCRS
)
{
const
int
loc
=
Utility
::
SearchSorted
(
_ja
,
j
,
_ia
[
i
],
_ia
[
i
+
1
]);
if
(
loc
>=
0
)
{
_val
[
loc
]
=
v
;
return
;
}
}
_tri
.
push_back
(
TRI_COORD
<
T
>
(
i
,
j
,
v
,
false
));
}
//-----------------------------------------------------------------------------
// Adds (i,j) to value
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
add
(
INDEX
i
,
INDEX
j
,
T
v
)
{
MICK
(
i
,
j
);
// Matrix Index ChecKing
if
(
i
<
_nRowsCRS
)
{
const
int
loc
=
Utility
::
SearchSorted
(
_ja
,
j
,
_ia
[
i
],
_ia
[
i
+
1
]);
if
(
loc
>=
0
)
{
_val
[
loc
]
+=
v
;
return
;
}
}
_tri
.
push_back
(
TRI_COORD
<
T
>
(
i
,
j
,
v
,
true
));
}
//-----------------------------------------------------------------------------
// returns a triplet value of the ith nonzero
//-----------------------------------------------------------------------------
template
<
typename
T
>
TRIPLET
<
T
>
SparseMatrix
<
T
>::
get_triplet
(
INDEX
i
)
const
{
compress
(
*
this
);
if
(
i
>=
_ia
[
_nRowsCRS
])
{
gerror
(
"ERROR: tried indexing triplet of sparse matrix beyond range"
);
}
INDEX
row
(
std
::
lower_bound
(
_ia
,
_ia
+
_nRowsCRS
,
i
)
-
_ia
);
row
-=
_ia
[
row
]
!=
i
;
return
TRIPLET
<
T
>
(
row
,
_ja
[
i
],
_val
[
i
]);
}
//-----------------------------------------------------------------------------
// full reset - completely wipes out all SparseMatrix data, zero is ignored
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
reset
(
INDEX
rows
,
INDEX
cols
,
bool
zero
)
{
_delete
();
_nRows
=
rows
;
_nCols
=
cols
;
}
//-----------------------------------------------------------------------------
// resize - changes the _nRows and _nCols without changing anything else
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
resize
(
INDEX
rows
,
INDEX
cols
,
bool
copy
)
{
GCHK
(
_nRowsCRS
>
rows
,
"SparseMatrix::resize CRS rows exceed specified rows"
);
_nRows
=
rows
;
_nCols
=
cols
;
// a check on this would be expensive
}
//-----------------------------------------------------------------------------
// get sparsity from DenseMatrix, if TOL < 0, then only zero values are added
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
reset
(
const
DenseMatrix
<
T
>&
D
,
double
TOL
)
{
_delete
();
// clears all values
// if TOL is specified then TOL = TOL^2 * max(abs(D))^2
if
(
TOL
>
0.0
)
{
TOL
*=
D
.
maxabs
();
TOL
*=
TOL
;
}
_nRows
=
D
.
nRows
();
_nCols
=
D
.
nCols
();
for
(
INDEX
i
=
0
;
i
<
D
.
nRows
();
i
++
)
for
(
INDEX
j
=
0
;
j
<
D
.
nCols
();
j
++
)
if
(
D
(
i
,
j
)
*
D
(
i
,
j
)
>=
TOL
)
// if TOL wasn't specified then TOL < 0
set
(
i
,
j
,
D
(
i
,
j
));
compress
();
}
//-----------------------------------------------------------------------------
// copy - dangerous: ignores rows & columns
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
copy
(
const
T
*
ptr
,
INDEX
rows
,
INDEX
cols
)
{
cout
<<
"SparseMatrix<T>::copy() has no effect.
\n
"
;
}
//-----------------------------------------------------------------------------
// dense_copy - copy to dense matrix
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
dense_copy
(
DenseMatrix
<
T
>
&
D
)
const
{
SparseMatrix
<
T
>::
compress
(
*
this
);
D
.
reset
(
nRows
(),
nCols
());
for
(
INDEX
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
for
(
INDEX
j
=
_ia
[
i
];
j
<
_ia
[
i
+
1
];
j
++
)
D
(
i
,
_ja
[
j
])
=
_val
[
j
];
}
template
<
typename
T
>
DenseMatrix
<
T
>
SparseMatrix
<
T
>::
dense_copy
(
void
)
const
{
DenseMatrix
<
T
>
D
;
dense_copy
(
D
);
return
D
;
}
//-----------------------------------------------------------------------------
// returns true if the matrix has no non-zero elements
//-----------------------------------------------------------------------------
template
<
typename
T
>
bool
SparseMatrix
<
T
>::
empty
()
const
{
return
_size
==
0
&&
_tri
.
empty
();
}
//-----------------------------------------------------------------------------
// returns the number of rows specified by the user
//-----------------------------------------------------------------------------
template
<
typename
T
>
inline
INDEX
SparseMatrix
<
T
>::
nRows
()
const
{
return
_nRows
;
}
//-----------------------------------------------------------------------------
// returns the number of columns specified by the user
//-----------------------------------------------------------------------------
template
<
typename
T
>
inline
INDEX
SparseMatrix
<
T
>::
nCols
()
const
{
return
_nCols
;
}
//-----------------------------------------------------------------------------
// returns the number of non-zeros in the matrix
//-----------------------------------------------------------------------------
template
<
typename
T
>
INDEX
SparseMatrix
<
T
>::
size
()
const
{
compress
(
*
this
);
return
_size
;
}
//-----------------------------------------------------------------------------
// returns the number of nonzero elements in a row
//-----------------------------------------------------------------------------
template
<
typename
T
>
INDEX
SparseMatrix
<
T
>::
RowSize
(
INDEX
r
)
const
{
compress
(
*
this
);
GCHK
(
r
>=
_nRows
,
"Rowsize: invalid row"
);
if
(
r
>=
_nRowsCRS
)
return
0
;
return
_ia
[
r
+
1
]
-
_ia
[
r
];
}
//-----------------------------------------------------------------------------
// returns a pointer to the data, causes a compress
//-----------------------------------------------------------------------------
template
<
typename
T
>
T
*
SparseMatrix
<
T
>::
get_ptr
()
const
{
compress
(
*
this
);
return
_val
;
}
//-----------------------------------------------------------------------------
// returns true if (i,j) falls in the user specified range
//-----------------------------------------------------------------------------
template
<
typename
T
>
bool
SparseMatrix
<
T
>::
in_range
(
INDEX
i
,
INDEX
j
)
const
{
return
i
<
nRows
()
&&
j
<
nCols
();
}
//-----------------------------------------------------------------------------
// assigns this sparsematrix from another one - full memory copy
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>&
SparseMatrix
<
T
>::
operator
=
(
const
SparseMatrix
<
T
>
&
C
)
{
_delete
();
_copy
(
C
);
return
*
this
;
}
//-----------------------------------------------------------------------------
// assigns this sparsematrix from another one - full memory copy
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>&
SparseMatrix
<
T
>::
operator
=
(
const
T
v
)
{
// set_all_elements only changes _data, so we need a compress
compress
(
*
this
);
return
set_all_elements_to
(
v
);
}
//-----------------------------------------------------------------------------
// scales this sparse matrix by a constant
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>&
SparseMatrix
<
T
>::
operator
*=
(
const
T
&
a
)
{
compress
(
*
this
);
for
(
INDEX
i
=
0
;
i
<
size
();
i
++
)
_val
[
i
]
*=
a
;
return
*
this
;
}
//-----------------------------------------------------------------------------
// Multiplies this Sparsematrix tranposed times a vector
//-----------------------------------------------------------------------------
template
<
typename
T
>
DenseVector
<
T
>
SparseMatrix
<
T
>::
transMat
(
const
Vector
<
T
>
&
x
)
const
{
DenseVector
<
T
>
y
(
nCols
(),
true
);
GCK
(
*
this
,
x
,
nRows
()
!=
x
.
size
(),
"operator *: Sparse matrix incompatible with Vector."
)
INDEX
i
,
ij
;
for
(
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
for
(
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
y
(
_ja
[
ij
])
+=
_val
[
ij
]
*
x
(
i
);
return
y
;
}
//-----------------------------------------------------------------------------
// Return matrix transpose
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>
SparseMatrix
<
T
>::
transpose
()
const
{
compress
(
*
this
);
SparseMatrix
<
T
>
At
(
nCols
(),
nRows
());
for
(
INDEX
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
for
(
INDEX
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
At
.
set
(
_ja
[
ij
],
i
,
_val
[
ij
]);
compress
(
At
);
return
At
;
}
//-----------------------------------------------------------------------------
// Matrix Transpose/DenseMatrix multiply
//-----------------------------------------------------------------------------
template
<
typename
T
>
DenseMatrix
<
T
>
SparseMatrix
<
T
>::
transMat
(
const
DenseMatrix
<
T
>
&
D
)
const
{
compress
(
*
this
);
GCK
(
*
this
,
D
,
nRows
()
!=
D
.
nRows
(),
"transMat: Sparse matrix incompatible with DenseMatrix."
)
DenseMatrix
<
T
>
C
(
nCols
(),
D
.
nCols
(),
true
);
// initialized to zero
INDEX
j
,
k
,
ki
;
for
(
k
=
0
;
k
<
_nRowsCRS
;
k
++
)
for
(
ki
=
_ia
[
k
];
ki
<
_ia
[
k
+
1
];
ki
++
)
for
(
j
=
0
;
j
<
D
.
nCols
();
j
++
)
C
(
_ja
[
ki
],
j
)
+=
_val
[
ki
]
*
D
(
k
,
j
);
// C(i,j) = S(k,i) * D(k, j)
return
C
;
}
//-----------------------------------------------------------------------------
// Matrix Transpose/SparseMatrix multiply - IS THIS REALLY NEEDED??
//-----------------------------------------------------------------------------
template
<
typename
T
>
DenseMatrix
<
T
>
SparseMatrix
<
T
>::
transMat
(
const
SparseMatrix
<
T
>
&
D
)
const
{
compress
(
*
this
);
GCK
(
*
this
,
D
,
nRows
()
!=
D
.
nRows
(),
"transMat: Sparse matrix incompatible with DenseMatrix."
)
DenseMatrix
<
T
>
C
(
nCols
(),
D
.
nCols
(),
true
);
// initialized to zero
INDEX
k
,
ki
,
kj
;
for
(
k
=
0
;
k
<
_nRowsCRS
;
k
++
)
for
(
kj
=
D
.
_ia
[
k
];
kj
<
D
.
_ia
[
k
+
1
];
kj
++
)
for
(
ki
=
_ia
[
k
];
ki
<
_ia
[
k
+
1
];
ki
++
)
C
(
_ja
[
ki
],
D
.
_ja
[
kj
])
+=
_val
[
ki
]
*
D
.
_val
[
kj
];
// C(i,j) = S(k,i)*D(k,j)
return
C
;
}
//-----------------------------------------------------------------------------
// multiplies each row by the corresponding element in Vector scale
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>&
SparseMatrix
<
T
>::
row_scale
(
const
Vector
<
T
>
&
v
)
{
compress
(
*
this
);
INDEX
i
,
ij
;
GCK
(
*
this
,
v
,
v
.
size
()
!=
nRows
(),
"Incompatible Vector length in row_scale."
);
for
(
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
for
(
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
_val
[
ij
]
*=
v
[
i
];
return
*
this
;
}
//-----------------------------------------------------------------------------
// multiples each column by the corresponding element in Vector scale
//-----------------------------------------------------------------------------
template
<
typename
T
>
SparseMatrix
<
T
>&
SparseMatrix
<
T
>::
col_scale
(
const
Vector
<
T
>
&
v
)
{
compress
(
*
this
);
INDEX
i
,
ij
;
GCK
(
*
this
,
v
,
v
.
size
()
!=
nCols
(),
"Incompatible Vector length in col_scale."
);
for
(
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
for
(
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
_val
[
ij
]
*=
v
[
_ja
[
ij
]];
return
*
this
;
}
//-----------------------------------------------------------------------------
// Returns a vector of the sums of each column
//-----------------------------------------------------------------------------
template
<
typename
T
>
DenseVector
<
T
>
SparseMatrix
<
T
>::
col_sum
()
const
{
compress
(
*
this
);
INDEX
i
,
ij
;
GCHK
(
!
nRows
(),
"SparseMatrix::Matrix not initialized in col_sum."
)
DenseVector
<
T
>
csum
(
nCols
());
for
(
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
for
(
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
csum
(
_ja
[
ij
])
+=
_val
[
ij
];
return
(
csum
);
}
//-----------------------------------------------------------------------------
// Returns a vector with the number of nonzeros in each column
//-----------------------------------------------------------------------------
template
<
typename
T
>
DenseVector
<
INDEX
>
SparseMatrix
<
T
>::
column_count
()
const
{
compress
(
*
this
);
INDEX
i
,
j
;
DenseVector
<
INDEX
>
counts
(
nCols
());
for
(
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
for
(
j
=
_ia
[
i
];
j
<
_ia
[
i
+
1
];
j
++
)
counts
(
_ja
[
j
])
++
;
return
(
counts
);
}
//-----------------------------------------------------------------------------
// Writes a the nonzeros of a row to a vector
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
get_row
(
INDEX
i
,
DenseVector
<
T
>&
row
,
DenseVector
<
INDEX
>&
indx
)
const
{
GCHK
(
i
>=
nRows
(),
"get_row() - invalid row number"
);
row
.
resize
(
RowSize
(
i
));
indx
.
resize
(
row
.
size
());
INDEX
idx
=
0
,
ij
;
for
(
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
{
row
(
idx
)
=
_val
[
ij
];
indx
(
idx
++
)
=
_ja
[
ij
];
}
}
//-----------------------------------------------------------------------------
// Computes the product of N'DN
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
WeightedLeastSquares
(
const
SparseMatrix
<
T
>
&
N
,
const
DiagonalMatrix
<
T
>
&
D
)
{
compress
(
N
);
GCK
(
N
,
D
,
N
.
nRows
()
!=
D
.
nRows
(),
"SparseMatrix::WeightedLeastSquares()"
);
INDEX
k
,
ki
,
kj
;
resize
(
N
.
nCols
(),
N
.
nCols
());
// set size of this matrix
for
(
k
=
0
;
k
<
_size
;
k
++
)
_val
[
k
]
=
0.0
;
// compute R(i,j) = N(k,i) D(k,q) N(i,j) = N(k,i)*D(k,k)*N(k,j) (sum on k)
for
(
k
=
0
;
k
<
N
.
_nRowsCRS
;
k
++
)
for
(
ki
=
N
.
_ia
[
k
];
ki
<
N
.
_ia
[
k
+
1
];
ki
++
)
for
(
kj
=
N
.
_ia
[
k
];
kj
<
N
.
_ia
[
k
+
1
];
kj
++
)
add
(
N
.
_ja
[
ki
],
N
.
_ja
[
kj
],
D
[
k
]
*
N
[
kj
]
*
N
[
ki
]);
compress
();
}
//-----------------------------------------------------------------------------
// Return a diagonal matrix containing the diagonal entries of this matrix
//-----------------------------------------------------------------------------
template
<
typename
T
>
DiagonalMatrix
<
T
>
SparseMatrix
<
T
>::
get_diag
()
const
{
compress
(
*
this
);
DiagonalMatrix
<
T
>
D
(
nRows
(),
true
);
// initialized to zero
INDEX
i
,
ij
;
for
(
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
{
for
(
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
{
if
(
_ja
[
ij
]
>=
i
)
// have we reached or passed the diagonal?
{
if
(
_ja
[
ij
]
==
i
)
D
[
i
]
=
_val
[
ij
];
// this this the diagonal?
break
;
// D[i] is already zero if there is no diagonal
}
}
}
return
D
;
}
//-----------------------------------------------------------------------------
// output function - builds a string with each nonzero triplet value
//-----------------------------------------------------------------------------
template
<
typename
T
>
string
SparseMatrix
<
T
>::
tostring
()
const
{
compress
(
*
this
);
string
out
;
INDEX
i
,
ij
;
for
(
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
{
for
(
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
{
if
(
ij
)
out
+=
"
\n
"
;
// append newline if not first nonzero
out
+=
"("
+
ATC_STRING
::
tostring
(
i
)
+
", "
;
// append "(i,"
out
+=
ATC_STRING
::
tostring
(
_ja
[
ij
])
+
") = "
;
// append "j) = "
out
+=
ATC_STRING
::
tostring
(
_val
[
ij
]);
// append "value"
}
}
return
out
;
// return the completed string
}
//-----------------------------------------------------------------------------
// returns the maximum value in the row
//-----------------------------------------------------------------------------
template
<
typename
T
>
T
SparseMatrix
<
T
>::
get_row_max
(
INDEX
row
)
const
{
compress
(
*
this
);
if
(
!
RowSize
(
row
))
return
(
T
)
0
;
// if there are no nonzeros in the row
INDEX
ij
;
T
max
=
_val
[
_ia
[
row
]];
for
(
ij
=
_ia
[
row
]
+
1
;
ij
<
_ia
[
row
+
1
];
ij
++
)
max
=
std
::
max
(
max
,
_val
[
ij
]);
return
max
;
}
//-----------------------------------------------------------------------------
// returns the minimum value in the row
//-----------------------------------------------------------------------------
template
<
typename
T
>
T
SparseMatrix
<
T
>::
get_row_min
(
INDEX
row
)
const
{
compress
(
*
this
);
if
(
!
RowSize
(
row
))
return
(
T
)
0
;
// if there are no nonzeros in the row
INDEX
ij
;
T
min
=
_val
[
_ia
[
row
]];
for
(
ij
=
_ia
[
row
]
+
1
;
ij
<
_ia
[
row
+
1
];
ij
++
)
min
=
std
::
min
(
min
,
_val
[
ij
]);
return
min
;
}
//-----------------------------------------------------------------------------
// prints a histogram of the values of a row to the screen
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
print_row_histogram
(
const
string
&
name
,
INDEX
nbins
)
const
{
compress
(
*
this
);
cout
<<
"Begin histogram "
<<
name
<<
"
\n
"
;
cout
<<
"# rows: "
<<
_nRows
<<
" columns: "
<<
_nCols
<<
" size: "
<<
_size
<<
"
\n
"
;
for
(
INDEX
i
=
0
;
i
<
_nRows
;
i
++
)
{
print_row_histogram
(
i
,
nbins
);
cout
<<
"
\n
"
;
}
cout
<<
"End histogram "
<<
name
<<
"
\n
"
;
}
//-----------------------------------------------------------------------------
// prints a histogram of the values of a row to the screen
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
print_row_histogram
(
INDEX
row
,
INDEX
nbins
)
const
{
compress
(
*
this
);
if
(
!
nbins
)
nbins
++
;
vector
<
INDEX
>
counts
(
nbins
,
0
);
const
T
min
=
get_row_min
(
row
);
const
T
max
=
get_row_max
(
row
);
const
T
range
=
max
-
min
;
const
double
bin_size
=
range
/
double
(
nbins
);
if
(
range
<=
0.0
)
counts
[
nbins
-
1
]
=
RowSize
(
row
);
else
{
for
(
INDEX
ij
=
_ia
[
row
];
ij
<
_ia
[
row
+
1
];
ij
++
)
{
INDEX
bin
=
INDEX
((
_val
[
ij
]
-
min
)
/
bin_size
);
counts
[
bin
-
(
bin
==
nbins
)]
++
;
}
}
cout
<<
showbase
<<
scientific
;
cout
<<
"# Histogram: row "
<<
row
<<
" min "
<<
min
<<
" max "
<<
max
<<
" cnt "
<<
RowSize
(
row
)
<<
"
\n
"
;
T
bin_start
=
min
;
for
(
INDEX
i
=
0
;
i
<
nbins
;
i
++
)
{
cout
<<
"("
<<
bin_start
<<
","
;
bin_start
+=
bin_size
;
cout
<<
bin_start
<<
") "
<<
counts
[
i
]
<<
"
\n
"
;
}
}
//-----------------------------------------------------------------------------
// Outputs a string to a sparse Matlab type
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
matlab
(
ostream
&
o
,
const
string
&
s
)
const
{
compress
(
*
this
);
INDEX
i
,
ij
;
o
<<
s
<<
" = sparse("
<<
nRows
()
<<
","
<<
nCols
()
<<
");
\n
"
;
o
<<
showbase
<<
scientific
;
for
(
i
=
0
;
i
<
_nRowsCRS
;
i
++
)
for
(
ij
=
_ia
[
i
];
ij
<
_ia
[
i
+
1
];
ij
++
)
o
<<
s
<<
"("
<<
i
+
1
<<
","
<<
_ja
[
ij
]
+
1
<<
")="
<<
_val
[
ij
]
<<
";
\n
"
;
}
//-----------------------------------------------------------------------------
// Writes the matrix to a binary file (after a compress).
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
binary_write
(
std
::
fstream
&
f
)
const
{
compress
(
*
this
);
f
.
write
((
char
*
)
&
_size
,
sizeof
(
INDEX
));
// writes number of nonzeros
f
.
write
((
char
*
)
&
_nRowsCRS
,
sizeof
(
INDEX
));
// writes number of rows in crs
f
.
write
((
char
*
)
&
_nRows
,
sizeof
(
INDEX
));
// write matrix rows
f
.
write
((
char
*
)
&
_nCols
,
sizeof
(
INDEX
));
// write number of columns
if
(
!
_size
)
return
;
f
.
write
((
char
*
)
_val
,
sizeof
(
T
)
*
_size
);
f
.
write
((
char
*
)
_ja
,
sizeof
(
INDEX
)
*
_size
);
f
.
write
((
char
*
)
_ia
,
sizeof
(
INDEX
)
*
(
_nRowsCRS
+
1
));
}
//-----------------------------------------------------------------------------
// Reads a SparseMatrix from a binary file. (wipes out any original data)
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
binary_read
(
std
::
fstream
&
f
)
{
_delete
();
f
.
read
((
char
*
)
&
_size
,
sizeof
(
INDEX
));
f
.
read
((
char
*
)
&
_nRowsCRS
,
sizeof
(
INDEX
));
f
.
read
((
char
*
)
&
_nRows
,
sizeof
(
INDEX
));
f
.
read
((
char
*
)
&
_nCols
,
sizeof
(
INDEX
));
if
(
!
_size
)
return
;
_create
(
_size
,
_nRowsCRS
);
f
.
read
((
char
*
)
_val
,
sizeof
(
T
)
*
_size
);
f
.
read
((
char
*
)
_ja
,
sizeof
(
INDEX
)
*
_size
);
f
.
read
((
char
*
)
_ia
,
sizeof
(
INDEX
)
*
(
_nRowsCRS
+
1
));
}
//-----------------------------------------------------------------------------
// Writes the sparse matrix to a file in a binary format
//-----------------------------------------------------------------------------
template
<
typename
T
>
void
SparseMatrix
<
T
>::
write_restart
(
FILE
*
f
)
const
{
compress
(
*
this
);
fwrite
(
&
_size
,
sizeof
(
INDEX
),
1
,
f
);
// write number of nonzeros
fwrite
(
&
_nRowsCRS
,
sizeof
(
INDEX
),
1
,
f
);
// write number of rows
fwrite
(
&
_nRows
,
sizeof
(
INDEX
),
1
,
f
);
// write number of columns
fwrite
(
&
_nCols
,
sizeof
(
INDEX
),
1
,
f
);
// write number of columns
if
(
!
_size
)
return
;
fwrite
(
_val
,
sizeof
(
T
),
_size
,
f
);
fwrite
(
_ja
,
sizeof
(
T
),
_size
,
f
);
fwrite
(
_ia
,
sizeof
(
INDEX
),
_nRowsCRS
+
1
,
f
);
}
#endif
Event Timeline
Log In to Comment