1 //=================================================================================================
2 /*!
3 // \file blaze/math/HermitianMatrix.h
4 // \brief Header file for the complete HermitianMatrix implementation
5 //
6 // Copyright (C) 2012-2020 Klaus Iglberger - All Rights Reserved
7 //
8 // This file is part of the Blaze library. You can redistribute it and/or modify it under
9 // the terms of the New (Revised) BSD License. Redistribution and use in source and binary
10 // forms, with or without modification, are permitted provided that the following conditions
11 // are met:
12 //
13 // 1. Redistributions of source code must retain the above copyright notice, this list of
14 // conditions and the following disclaimer.
15 // 2. Redistributions in binary form must reproduce the above copyright notice, this list
16 // of conditions and the following disclaimer in the documentation and/or other materials
17 // provided with the distribution.
18 // 3. Neither the names of the Blaze development group nor the names of its contributors
19 // may be used to endorse or promote products derived from this software without specific
20 // prior written permission.
21 //
22 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
23 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
24 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
25 // SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
27 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
28 // BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 // ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
31 // DAMAGE.
32 */
33 //=================================================================================================
34
35 #ifndef _BLAZE_MATH_HERMITIANMATRIX_H_
36 #define _BLAZE_MATH_HERMITIANMATRIX_H_
37
38
39 //*************************************************************************************************
40 // Includes
41 //*************************************************************************************************
42
43 #include <cmath>
44 #include <vector>
45 #include <blaze/math/Aliases.h>
46 #include <blaze/math/adaptors/DiagonalMatrix.h>
47 #include <blaze/math/adaptors/HermitianMatrix.h>
48 #include <blaze/math/constraints/DenseMatrix.h>
49 #include <blaze/math/constraints/Resizable.h>
50 #include <blaze/math/constraints/SparseMatrix.h>
51 #include <blaze/math/DenseMatrix.h>
52 #include <blaze/math/Exception.h>
53 #include <blaze/math/shims/Real.h>
54 #include <blaze/math/SparseMatrix.h>
55 #include <blaze/math/typetraits/IsDenseMatrix.h>
56 #include <blaze/math/typetraits/UnderlyingBuiltin.h>
57 #include <blaze/util/Assert.h>
58 #include <blaze/util/IntegralConstant.h>
59 #include <blaze/util/Random.h>
60 #include <blaze/util/Types.h>
61
62
63 namespace blaze {
64
65 //=================================================================================================
66 //
67 // RAND SPECIALIZATION
68 //
69 //=================================================================================================
70
71 //*************************************************************************************************
72 /*! \cond BLAZE_INTERNAL */
73 /*!\brief Specialization of the Rand class template for HermitianMatrix.
74 // \ingroup random
75 //
76 // This specialization of the Rand class creates random instances of HermitianMatrix.
77 */
78 template< typename MT // Type of the adapted matrix
79 , bool SO // Storage order of the adapted matrix
80 , bool DF > // Density flag
81 class Rand< HermitianMatrix<MT,SO,DF> >
82 {
83 public:
84 //**********************************************************************************************
85 /*!\brief Generation of a random HermitianMatrix.
86 //
87 // \return The generated random matrix.
88 */
generate()89 inline const HermitianMatrix<MT,SO,DF> generate() const
90 {
91 BLAZE_CONSTRAINT_MUST_NOT_BE_RESIZABLE_TYPE( MT );
92
93 HermitianMatrix<MT,SO,DF> matrix;
94 randomize( matrix );
95 return matrix;
96 }
97 //**********************************************************************************************
98
99 //**********************************************************************************************
100 /*!\brief Generation of a random HermitianMatrix.
101 //
102 // \param n The number of rows and columns of the random matrix.
103 // \return The generated random matrix.
104 */
generate(size_t n)105 inline const HermitianMatrix<MT,SO,DF> generate( size_t n ) const
106 {
107 BLAZE_CONSTRAINT_MUST_BE_RESIZABLE_TYPE( MT );
108
109 HermitianMatrix<MT,SO,DF> matrix( n );
110 randomize( matrix );
111 return matrix;
112 }
113 //**********************************************************************************************
114
115 //**********************************************************************************************
116 /*!\brief Generation of a random HermitianMatrix.
117 //
118 // \param n The number of rows and columns of the random matrix.
119 // \param nonzeros The number of non-zero elements of the random matrix.
120 // \return The generated random matrix.
121 // \exception std::invalid_argument Invalid number of non-zero elements.
122 */
generate(size_t n,size_t nonzeros)123 inline const HermitianMatrix<MT,SO,DF> generate( size_t n, size_t nonzeros ) const
124 {
125 BLAZE_CONSTRAINT_MUST_BE_RESIZABLE_TYPE ( MT );
126 BLAZE_CONSTRAINT_MUST_BE_SPARSE_MATRIX_TYPE( MT );
127
128 if( nonzeros > n*n ) {
129 BLAZE_THROW_INVALID_ARGUMENT( "Invalid number of non-zero elements" );
130 }
131
132 HermitianMatrix<MT,SO,DF> matrix( n );
133 randomize( matrix, nonzeros );
134
135 return matrix;
136 }
137 //**********************************************************************************************
138
139 //**********************************************************************************************
140 /*!\brief Generation of a random HermitianMatrix.
141 //
142 // \param min The smallest possible value for a matrix element.
143 // \param max The largest possible value for a matrix element.
144 // \return The generated random matrix.
145 */
146 template< typename Arg > // Min/max argument type
generate(const Arg & min,const Arg & max)147 inline const HermitianMatrix<MT,SO,DF> generate( const Arg& min, const Arg& max ) const
148 {
149 BLAZE_CONSTRAINT_MUST_NOT_BE_RESIZABLE_TYPE( MT );
150
151 HermitianMatrix<MT,SO,DF> matrix;
152 randomize( matrix, min, max );
153 return matrix;
154 }
155 //**********************************************************************************************
156
157 //**********************************************************************************************
158 /*!\brief Generation of a random HermitianMatrix.
159 //
160 // \param n The number of rows and columns of the random matrix.
161 // \param min The smallest possible value for a matrix element.
162 // \param max The largest possible value for a matrix element.
163 // \return The generated random matrix.
164 */
165 template< typename Arg > // Min/max argument type
generate(size_t n,const Arg & min,const Arg & max)166 inline const HermitianMatrix<MT,SO,DF> generate( size_t n, const Arg& min, const Arg& max ) const
167 {
168 BLAZE_CONSTRAINT_MUST_BE_RESIZABLE_TYPE( MT );
169
170 HermitianMatrix<MT,SO,DF> matrix( n );
171 randomize( matrix, min, max );
172 return matrix;
173 }
174 //**********************************************************************************************
175
176 //**********************************************************************************************
177 /*!\brief Generation of a random HermitianMatrix.
178 //
179 // \param n The number of rows and columns of the random matrix.
180 // \param nonzeros The number of non-zero elements of the random matrix.
181 // \param min The smallest possible value for a matrix element.
182 // \param max The largest possible value for a matrix element.
183 // \return The generated random matrix.
184 // \exception std::invalid_argument Invalid number of non-zero elements.
185 */
186 template< typename Arg > // Min/max argument type
generate(size_t n,size_t nonzeros,const Arg & min,const Arg & max)187 inline const HermitianMatrix<MT,SO,DF> generate( size_t n, size_t nonzeros,
188 const Arg& min, const Arg& max ) const
189 {
190 BLAZE_CONSTRAINT_MUST_BE_RESIZABLE_TYPE ( MT );
191 BLAZE_CONSTRAINT_MUST_BE_SPARSE_MATRIX_TYPE( MT );
192
193 if( nonzeros > n*n ) {
194 BLAZE_THROW_INVALID_ARGUMENT( "Invalid number of non-zero elements" );
195 }
196
197 HermitianMatrix<MT,SO,DF> matrix( n );
198 randomize( matrix, nonzeros, min, max );
199
200 return matrix;
201 }
202 //**********************************************************************************************
203
204 //**********************************************************************************************
205 /*!\brief Randomization of a HermitianMatrix.
206 //
207 // \param matrix The matrix to be randomized.
208 // \return void
209 */
randomize(HermitianMatrix<MT,SO,DF> & matrix)210 inline void randomize( HermitianMatrix<MT,SO,DF>& matrix ) const
211 {
212 randomize( matrix, typename IsDenseMatrix<MT>::Type() );
213 }
214 //**********************************************************************************************
215
216 //**********************************************************************************************
217 /*!\brief Randomization of a sparse HermitianMatrix.
218 //
219 // \param matrix The matrix to be randomized.
220 // \param nonzeros The number of non-zero elements of the random matrix.
221 // \return void
222 // \exception std::invalid_argument Invalid number of non-zero elements.
223 */
randomize(HermitianMatrix<MT,SO,DF> & matrix,size_t nonzeros)224 inline void randomize( HermitianMatrix<MT,SO,DF>& matrix, size_t nonzeros ) const
225 {
226 BLAZE_CONSTRAINT_MUST_BE_SPARSE_MATRIX_TYPE( MT );
227
228 using ET = ElementType_t<MT>;
229 using BT = UnderlyingBuiltin_t<ET>;
230
231 const size_t n( matrix.rows() );
232
233 if( nonzeros > n*n ) {
234 BLAZE_THROW_INVALID_ARGUMENT( "Invalid number of non-zero elements" );
235 }
236
237 if( n == 0UL ) return;
238
239 matrix.reset();
240 matrix.reserve( nonzeros );
241
242 while( matrix.nonZeros() < nonzeros )
243 {
244 const size_t row ( rand<size_t>( 0UL, n-1UL ) );
245 const size_t column( rand<size_t>( 0UL, n-1UL ) );
246
247 if( row == column )
248 matrix(row,column) = rand<BT>();
249 else
250 matrix(row,column) = rand<ET>();
251 }
252 }
253 //**********************************************************************************************
254
255 //**********************************************************************************************
256 /*!\brief Randomization of a HermitianMatrix.
257 //
258 // \param matrix The matrix to be randomized.
259 // \param min The smallest possible value for a matrix element.
260 // \param max The largest possible value for a matrix element.
261 // \return void
262 */
263 template< typename Arg > // Min/max argument type
randomize(HermitianMatrix<MT,SO,DF> & matrix,const Arg & min,const Arg & max)264 inline void randomize( HermitianMatrix<MT,SO,DF>& matrix, const Arg& min, const Arg& max ) const
265 {
266 randomize( matrix, min, max, typename IsDenseMatrix<MT>::Type() );
267 }
268 //**********************************************************************************************
269
270 //**********************************************************************************************
271 /*!\brief Randomization of a sparse HermitianMatrix.
272 //
273 // \param matrix The matrix to be randomized.
274 // \param nonzeros The number of non-zero elements of the random matrix.
275 // \param min The smallest possible value for a matrix element.
276 // \param max The largest possible value for a matrix element.
277 // \return void
278 // \exception std::invalid_argument Invalid number of non-zero elements.
279 */
280 template< typename Arg > // Min/max argument type
randomize(HermitianMatrix<MT,SO,DF> & matrix,size_t nonzeros,const Arg & min,const Arg & max)281 inline void randomize( HermitianMatrix<MT,SO,DF>& matrix,
282 size_t nonzeros, const Arg& min, const Arg& max ) const
283 {
284 BLAZE_CONSTRAINT_MUST_BE_SPARSE_MATRIX_TYPE( MT );
285
286 using ET = ElementType_t<MT>;
287 using BT = UnderlyingBuiltin_t<ET>;
288
289 const size_t n( matrix.rows() );
290
291 if( nonzeros > n*n ) {
292 BLAZE_THROW_INVALID_ARGUMENT( "Invalid number of non-zero elements" );
293 }
294
295 if( n == 0UL ) return;
296
297 std::vector<size_t> dist( n );
298 std::vector<bool> structure( n*n );
299 size_t nz( 0UL );
300
301 while( nz < nonzeros )
302 {
303 const size_t row = rand<size_t>( 0UL, n-1UL );
304 const size_t col = rand<size_t>( 0UL, n-1UL );
305
306 if( structure[row*n+col] ) continue;
307
308 ++dist[row];
309 structure[row*n+col] = true;
310 ++nz;
311
312 if( row != col ) {
313 ++dist[col];
314 structure[col*n+row] = true;
315 ++nz;
316 }
317 }
318
319 matrix.reset();
320 matrix.reserve( nz );
321
322 for( size_t i=0UL; i<n; ++i ) {
323 matrix.reserve( i, dist[i] );
324 }
325
326 for( size_t i=0UL; i<n; ++i ) {
327 for( size_t j=i; j<n; ++j ) {
328 if( structure[i*n+j] ) {
329 if( i == j )
330 matrix.append( i, j, rand<BT>( real( min ), real( max ) ) );
331 else
332 matrix.append( i, j, rand<ET>( min, max ) );
333 }
334 }
335 }
336 }
337 //**********************************************************************************************
338
339 private:
340 //**********************************************************************************************
341 /*!\brief Randomization of a dense HermitianMatrix.
342 //
343 // \param matrix The matrix to be randomized.
344 // \return void
345 */
randomize(HermitianMatrix<MT,SO,DF> & matrix,TrueType)346 inline void randomize( HermitianMatrix<MT,SO,DF>& matrix, TrueType ) const
347 {
348 BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( MT );
349
350 using ET = ElementType_t<MT>;
351 using BT = UnderlyingBuiltin_t<ET>;
352
353 const size_t n( matrix.rows() );
354
355 for( size_t i=0UL; i<n; ++i ) {
356 for( size_t j=0UL; j<i; ++j ) {
357 matrix(i,j) = rand<ET>();
358 }
359 matrix(i,i) = rand<BT>();
360 }
361 }
362 //**********************************************************************************************
363
364 //**********************************************************************************************
365 /*!\brief Randomization of a sparse HermitianMatrix.
366 //
367 // \param matrix The matrix to be randomized.
368 // \return void
369 */
randomize(HermitianMatrix<MT,SO,DF> & matrix,FalseType)370 inline void randomize( HermitianMatrix<MT,SO,DF>& matrix, FalseType ) const
371 {
372 BLAZE_CONSTRAINT_MUST_BE_SPARSE_MATRIX_TYPE( MT );
373
374 const size_t n( matrix.rows() );
375
376 if( n == 0UL ) return;
377
378 const size_t nonzeros( rand<size_t>( 1UL, std::ceil( 0.5*n*n ) ) );
379
380 randomize( matrix, nonzeros );
381 }
382 //**********************************************************************************************
383
384 //**********************************************************************************************
385 /*!\brief Randomization of a dense HermitianMatrix.
386 //
387 // \param matrix The matrix to be randomized.
388 // \param min The smallest possible value for a matrix element.
389 // \param max The largest possible value for a matrix element.
390 // \return void
391 */
392 template< typename Arg > // Min/max argument type
randomize(HermitianMatrix<MT,SO,DF> & matrix,const Arg & min,const Arg & max,TrueType)393 inline void randomize( HermitianMatrix<MT,SO,DF>& matrix,
394 const Arg& min, const Arg& max, TrueType ) const
395 {
396 BLAZE_CONSTRAINT_MUST_BE_DENSE_MATRIX_TYPE( MT );
397
398 using ET = ElementType_t<MT>;
399 using BT = UnderlyingBuiltin_t<ET>;
400
401 const size_t n( matrix.rows() );
402
403 for( size_t i=0UL; i<n; ++i ) {
404 for( size_t j=0UL; j<i; ++j ) {
405 matrix(i,j) = rand<ET>( min, max );
406 }
407 matrix(i,i) = rand<BT>( real( min ), real( max ) );
408 }
409 }
410 //**********************************************************************************************
411
412 //**********************************************************************************************
413 /*!\brief Randomization of a sparse HermitianMatrix.
414 //
415 // \param matrix The matrix to be randomized.
416 // \param min The smallest possible value for a matrix element.
417 // \param max The largest possible value for a matrix element.
418 // \return void
419 */
420 template< typename Arg > // Min/max argument type
randomize(HermitianMatrix<MT,SO,DF> & matrix,const Arg & min,const Arg & max,FalseType)421 inline void randomize( HermitianMatrix<MT,SO,DF>& matrix,
422 const Arg& min, const Arg& max, FalseType ) const
423 {
424 BLAZE_CONSTRAINT_MUST_BE_SPARSE_MATRIX_TYPE( MT );
425
426 const size_t n( matrix.rows() );
427
428 if( n == 0UL ) return;
429
430 const size_t nonzeros( rand<size_t>( 1UL, std::ceil( 0.5*n*n ) ) );
431
432 randomize( matrix, nonzeros, min, max );
433 }
434 //**********************************************************************************************
435 };
436 /*! \endcond */
437 //*************************************************************************************************
438
439
440
441
442 //=================================================================================================
443 //
444 // MAKE FUNCTIONS
445 //
446 //=================================================================================================
447
448 //*************************************************************************************************
449 /*! \cond BLAZE_INTERNAL */
450 /*!\brief Setup of a random symmetric HermitianMatrix.
451 //
452 // \param matrix The matrix to be randomized.
453 // \return void
454 */
455 template< typename MT // Type of the adapted matrix
456 , bool SO // Storage order of the adapted matrix
457 , bool DF > // Density flag
makeSymmetric(HermitianMatrix<MT,SO,DF> & matrix)458 void makeSymmetric( HermitianMatrix<MT,SO,DF>& matrix )
459 {
460 using BT = UnderlyingBuiltin_t< ElementType_t<MT> >;
461
462 const size_t n( matrix.rows() );
463
464 for( size_t i=0UL; i<n; ++i ) {
465 for( size_t j=0UL; j<=i; ++j ) {
466 matrix(i,j) = rand<BT>();
467 }
468 }
469
470 BLAZE_INTERNAL_ASSERT( isSymmetric( matrix ), "Non-symmetric matrix detected" );
471 }
472 /*! \endcond */
473 //*************************************************************************************************
474
475
476 //*************************************************************************************************
477 /*! \cond BLAZE_INTERNAL */
478 /*!\brief Setup of a random symmetric HermitianMatrix.
479 //
480 // \param matrix The matrix to be randomized.
481 // \param min The smallest possible value for a matrix element.
482 // \param max The largest possible value for a matrix element.
483 // \return void
484 */
485 template< typename MT // Type of the adapted matrix
486 , bool SO // Storage order of the adapted matrix
487 , bool DF // Density flag
488 , typename Arg > // Min/max argument type
makeSymmetric(HermitianMatrix<MT,SO,DF> & matrix,const Arg & min,const Arg & max)489 void makeSymmetric( HermitianMatrix<MT,SO,DF>& matrix, const Arg& min, const Arg& max )
490 {
491 using BT = UnderlyingBuiltin_t< ElementType_t<MT> >;
492
493 const size_t n( matrix.rows() );
494
495 for( size_t i=0UL; i<n; ++i ) {
496 for( size_t j=0UL; j<=i; ++j ) {
497 matrix(i,j) = rand<BT>( real( min ), real( max ) );
498 }
499 }
500
501 BLAZE_INTERNAL_ASSERT( isSymmetric( matrix ), "Non-symmetric matrix detected" );
502 }
503 /*! \endcond */
504 //*************************************************************************************************
505
506
507 //*************************************************************************************************
508 /*! \cond BLAZE_INTERNAL */
509 /*!\brief Setup of a random Hermitian HermitianMatrix.
510 //
511 // \param matrix The matrix to be randomized.
512 // \return void
513 */
514 template< typename MT // Type of the adapted matrix
515 , bool SO // Storage order of the adapted matrix
516 , bool DF > // Density flag
makeHermitian(HermitianMatrix<MT,SO,DF> & matrix)517 void makeHermitian( HermitianMatrix<MT,SO,DF>& matrix )
518 {
519 using blaze::randomize;
520
521 randomize( matrix );
522 }
523 /*! \endcond */
524 //*************************************************************************************************
525
526
527 //*************************************************************************************************
528 /*! \cond BLAZE_INTERNAL */
529 /*!\brief Setup of a random Hermitian HermitianMatrix.
530 //
531 // \param matrix The matrix to be randomized.
532 // \param min The smallest possible value for a matrix element.
533 // \param max The largest possible value for a matrix element.
534 // \return void
535 */
536 template< typename MT // Type of the adapted matrix
537 , bool SO // Storage order of the adapted matrix
538 , bool DF // Density flag
539 , typename Arg > // Min/max argument type
makeHermitian(HermitianMatrix<MT,SO,DF> & matrix,const Arg & min,const Arg & max)540 void makeHermitian( HermitianMatrix<MT,SO,DF>& matrix, const Arg& min, const Arg& max )
541 {
542 using blaze::randomize;
543
544 randomize( matrix, min, max );
545 }
546 /*! \endcond */
547 //*************************************************************************************************
548
549
550 //*************************************************************************************************
551 /*! \cond BLAZE_INTERNAL */
552 /*!\brief Setup of a random (Hermitian) positive definite HermitianMatrix.
553 //
554 // \param matrix The matrix to be randomized.
555 // \return void
556 */
557 template< typename MT // Type of the adapted matrix
558 , bool SO // Storage order of the adapted matrix
559 , bool DF > // Density flag
makePositiveDefinite(HermitianMatrix<MT,SO,DF> & matrix)560 void makePositiveDefinite( HermitianMatrix<MT,SO,DF>& matrix )
561 {
562 using blaze::randomize;
563
564 using BT = UnderlyingBuiltin_t< ElementType_t<MT> >;
565
566 const size_t n( matrix.rows() );
567
568 randomize( matrix );
569 matrix *= matrix;
570
571 for( size_t i=0UL; i<n; ++i ) {
572 matrix(i,i) += BT(n);
573 }
574 }
575 /*! \endcond */
576 //*************************************************************************************************
577
578 } // namespace blaze
579
580 #endif
581