blob: 1f4a166a44d2cf8a6d740c77841cc81281320713 [file] [log] [blame]
Gael Guennebaud86ccd992008-11-05 13:47:55 +00001// This file is part of Eigen, a lightweight C++ template library
Benoit Jacob6347b1d2009-05-22 20:25:33 +02002// for linear algebra.
Gael Guennebaud86ccd992008-11-05 13:47:55 +00003//
Gael Guennebaud22c76092011-03-22 11:58:22 +01004// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
Gael Guennebaud86ccd992008-11-05 13:47:55 +00005//
Benoit Jacob69124cf2012-07-13 14:42:47 -04006// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
Gael Guennebaud86ccd992008-11-05 13:47:55 +00009
10#ifndef EIGEN_TESTSPARSE_H
Gael Guennebaud70df09b2011-10-24 09:31:33 +020011#define EIGEN_TESTSPARSE_H
Gael Guennebaud86ccd992008-11-05 13:47:55 +000012
Gael Guennebaud32124bc2011-01-27 17:36:58 +010013#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
14
Gael Guennebaud6a722602009-01-23 12:26:32 +000015#include "main.h"
16
Gael Guennebaud42e25782011-08-19 14:18:05 +020017#ifdef min
18#undef min
19#endif
20
21#ifdef max
22#undef max
23#endif
24
Gael Guennebaud4e8047c2019-02-20 13:59:34 +010025#include <unordered_map>
Gael Guennebaud6a722602009-01-23 12:26:32 +000026#define EIGEN_UNORDERED_MAP_SUPPORT
Gael Guennebaud4e8047c2019-02-20 13:59:34 +010027
Gael Guennebaud86ccd992008-11-05 13:47:55 +000028#include <Eigen/Cholesky>
29#include <Eigen/LU>
30#include <Eigen/Sparse>
31
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000032enum { ForceNonZeroDiag = 1, MakeLowerTriangular = 2, MakeUpperTriangular = 4, ForceRealDiag = 8 };
Gael Guennebaud86ccd992008-11-05 13:47:55 +000033
34/* Initializes both a sparse and dense matrix with same random values,
35 * and a ratio of \a density non zero entries.
36 * \param flags is a union of ForceNonZeroDiag, MakeLowerTriangular and MakeUpperTriangular
37 * allowing to control the shape of the matrix.
38 * \param zeroCoords and nonzeroCoords allows to get the coordinate lists of the non zero,
39 * and zero coefficients respectively.
40 */
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000041template <typename Scalar, int Opt1, int Opt2, typename StorageIndex>
42void initSparse(double density, Matrix<Scalar, Dynamic, Dynamic, Opt1>& refMat,
43 SparseMatrix<Scalar, Opt2, StorageIndex>& sparseMat, int flags = 0,
44 std::vector<Matrix<StorageIndex, 2, 1> >* zeroCoords = 0,
45 std::vector<Matrix<StorageIndex, 2, 1> >* nonzeroCoords = 0) {
46 enum { IsRowMajor = SparseMatrix<Scalar, Opt2, StorageIndex>::IsRowMajor };
Gael Guennebaud28293142009-05-04 14:25:12 +000047 sparseMat.setZero();
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000048 // sparseMat.reserve(int(refMat.rows()*refMat.cols()*density));
Erik Schultheisd271a7d2022-01-26 18:16:19 +000049 int nnz = static_cast<int>((1.5 * density) * static_cast<double>(IsRowMajor ? refMat.cols() : refMat.rows()));
50 sparseMat.reserve(VectorXi::Constant(IsRowMajor ? refMat.rows() : refMat.cols(), nnz));
Erik Schultheisb0fb5412021-11-18 18:33:31 +000051
52 Index insert_count = 0;
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000053 for (Index j = 0; j < sparseMat.outerSize(); j++) {
54 // sparseMat.startVec(j);
55 for (Index i = 0; i < sparseMat.innerSize(); i++) {
Gael Guennebaudb47ef142014-07-08 16:47:11 +020056 Index ai(i), aj(j);
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000057 if (IsRowMajor) std::swap(ai, aj);
58 Scalar v = (internal::random<double>(0, 1) < density) ? internal::random<Scalar>() : Scalar(0);
59 if ((flags & ForceNonZeroDiag) && (i == j)) {
Gael Guennebaudb986c142015-08-04 16:12:16 +020060 // FIXME: the following is too conservative
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000061 v = internal::random<Scalar>() * Scalar(3.);
62 v = v * v;
63 if (numext::real(v) > 0)
64 v += Scalar(5);
65 else
66 v -= Scalar(5);
Gael Guennebaud86ccd992008-11-05 13:47:55 +000067 }
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000068 if ((flags & MakeLowerTriangular) && aj > ai)
Gael Guennebaud86ccd992008-11-05 13:47:55 +000069 v = Scalar(0);
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000070 else if ((flags & MakeUpperTriangular) && aj < ai)
Gael Guennebaud86ccd992008-11-05 13:47:55 +000071 v = Scalar(0);
Gael Guennebaud8710bd22010-06-02 13:32:13 +020072
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000073 if ((flags & ForceRealDiag) && (i == j)) v = numext::real(v);
Gael Guennebaud8710bd22010-06-02 13:32:13 +020074
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000075 if (!numext::is_exactly_zero(v)) {
76 // sparseMat.insertBackByOuterInner(j,i) = v;
77 sparseMat.insertByOuterInner(j, i) = v;
Erik Schultheisb0fb5412021-11-18 18:33:31 +000078 ++insert_count;
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000079 if (nonzeroCoords) nonzeroCoords->push_back(Matrix<StorageIndex, 2, 1>(ai, aj));
80 } else if (zeroCoords) {
81 zeroCoords->push_back(Matrix<StorageIndex, 2, 1>(ai, aj));
Gael Guennebaud86ccd992008-11-05 13:47:55 +000082 }
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000083 refMat(ai, aj) = v;
Erik Schultheisb0fb5412021-11-18 18:33:31 +000084
85 // make sure we only insert as many as the sparse matrix supports
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000086 if (insert_count == NumTraits<StorageIndex>::highest()) return;
Gael Guennebaud86ccd992008-11-05 13:47:55 +000087 }
88 }
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000089 // sparseMat.finalize();
Gael Guennebaud86ccd992008-11-05 13:47:55 +000090}
91
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000092template <typename Scalar, int Options, typename Index>
93void initSparse(double density, Matrix<Scalar, Dynamic, 1>& refVec, SparseVector<Scalar, Options, Index>& sparseVec,
94 std::vector<int>* zeroCoords = 0, std::vector<int>* nonzeroCoords = 0) {
95 sparseVec.reserve(int(refVec.size() * density));
Gael Guennebaud709e9032009-01-07 17:01:57 +000096 sparseVec.setZero();
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +000097 for (int i = 0; i < refVec.size(); i++) {
98 Scalar v = (internal::random<double>(0, 1) < density) ? internal::random<Scalar>() : Scalar(0);
99 if (!numext::is_exactly_zero(v)) {
Gael Guennebaud28293142009-05-04 14:25:12 +0000100 sparseVec.insertBack(i) = v;
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +0000101 if (nonzeroCoords) nonzeroCoords->push_back(i);
102 } else if (zeroCoords)
103 zeroCoords->push_back(i);
Gael Guennebaud709e9032009-01-07 17:01:57 +0000104 refVec[i] = v;
105 }
106}
107
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +0000108template <typename Scalar, int Options, typename Index>
109void initSparse(double density, Matrix<Scalar, 1, Dynamic>& refVec, SparseVector<Scalar, Options, Index>& sparseVec,
110 std::vector<int>* zeroCoords = 0, std::vector<int>* nonzeroCoords = 0) {
111 sparseVec.reserve(int(refVec.size() * density));
Gael Guennebaud98ce4452013-03-06 11:58:22 +0100112 sparseVec.setZero();
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +0000113 for (int i = 0; i < refVec.size(); i++) {
114 Scalar v = (internal::random<double>(0, 1) < density) ? internal::random<Scalar>() : Scalar(0);
115 if (v != Scalar(0)) {
Gael Guennebaud98ce4452013-03-06 11:58:22 +0100116 sparseVec.insertBack(i) = v;
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +0000117 if (nonzeroCoords) nonzeroCoords->push_back(i);
118 } else if (zeroCoords)
119 zeroCoords->push_back(i);
Gael Guennebaud98ce4452013-03-06 11:58:22 +0100120 refVec[i] = v;
121 }
122}
123
Antonio Sánchez46e9cdb2023-12-05 21:22:55 +0000124#endif // EIGEN_TESTSPARSE_H