File: matrix_la_abstract.h

package info (click to toggle)
mldemos 0.5.1-3
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 32,224 kB
  • ctags: 46,525
  • sloc: cpp: 306,887; ansic: 167,718; ml: 126; sh: 109; makefile: 2
file content (840 lines) | stat: -rw-r--r-- 32,043 bytes parent folder | download | duplicates (4)
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
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
// Copyright (C) 2009  Davis E. King (davis@dlib.net)
// License: Boost Software License   See LICENSE.txt for the full license.
#undef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_
#ifdef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_ 

#include "matrix_abstract.h"
#include <complex>

namespace dlib
{

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
//                             Global linear algebra functions 
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    const matrix_exp::matrix_type inv (
        const matrix_exp& m
    );
    /*!
        requires
            - m is a square matrix
        ensures
            - returns the inverse of m 
              (Note that if m is singular or so close to being singular that there
              is a lot of numerical error then the returned matrix will be bogus.  
              You can check by seeing if m*inv(m) is an identity matrix)
    !*/

// ----------------------------------------------------------------------------------------

    const matrix pinv (
        const matrix_exp& m
    );
    /*!
        ensures
            - returns the Moore-Penrose pseudoinverse of m.
            - The returned matrix has m.nc() rows and m.nr() columns.
    !*/

// ----------------------------------------------------------------------------------------

    void svd (
        const matrix_exp& m,
        matrix<matrix_exp::type>& u,
        matrix<matrix_exp::type>& w,
        matrix<matrix_exp::type>& v
    );
    /*!
        ensures
            - computes the singular value decomposition of m
            - m == #u*#w*trans(#v)
            - trans(#u)*#u == identity matrix
            - trans(#v)*#v == identity matrix
            - diag(#w) == the singular values of the matrix m in no 
              particular order.  All non-diagonal elements of #w are
              set to 0.
            - #u.nr() == m.nr()
            - #u.nc() == m.nc()
            - #w.nr() == m.nc()
            - #w.nc() == m.nc()
            - #v.nr() == m.nc()
            - #v.nc() == m.nc()
            - if DLIB_USE_LAPACK is #defined then the xGESVD routine
              from LAPACK is used to compute the SVD.
    !*/

// ----------------------------------------------------------------------------------------

    long svd2 (
        bool withu, 
        bool withv, 
        const matrix_exp& m,
        matrix<matrix_exp::type>& u,
        matrix<matrix_exp::type>& w,
        matrix<matrix_exp::type>& v
    );
    /*!
        requires
            - m.nr() >= m.nc()
        ensures
            - computes the singular value decomposition of matrix m
            - m == subm(#u,get_rect(m))*diagm(#w)*trans(#v)
            - trans(#u)*#u == identity matrix
            - trans(#v)*#v == identity matrix
            - #w == the singular values of the matrix m in no 
              particular order.  
            - #u.nr() == m.nr()
            - #u.nc() == m.nr()
            - #w.nr() == m.nc()
            - #w.nc() == 1 
            - #v.nr() == m.nc()
            - #v.nc() == m.nc()
            - if (widthu == false) then
                - ignore the above regarding #u, it isn't computed and its
                  output state is undefined.
            - if (widthv == false) then
                - ignore the above regarding #v, it isn't computed and its
                  output state is undefined.
            - returns an error code of 0, if no errors and 'k' if we fail to
              converge at the 'kth' singular value.
            - if (DLIB_USE_LAPACK is #defined) then 
                - if (withu == withv) then
                    - the xGESDD routine from LAPACK is used to compute the SVD.
                - else
                    - the xGESVD routine from LAPACK is used to compute the SVD.
    !*/

// ----------------------------------------------------------------------------------------

    void svd3 (
        const matrix_exp& m,
        matrix<matrix_exp::type>& u,
        matrix<matrix_exp::type>& w,
        matrix<matrix_exp::type>& v
    );
    /*!
        ensures
            - computes the singular value decomposition of m
            - m == #u*diagm(#w)*trans(#v)
            - trans(#u)*#u == identity matrix
            - trans(#v)*#v == identity matrix
            - #w == the singular values of the matrix m in no 
              particular order.  
            - #u.nr() == m.nr()
            - #u.nc() == m.nc()
            - #w.nr() == m.nc()
            - #w.nc() == 1 
            - #v.nr() == m.nc()
            - #v.nc() == m.nc()
            - if DLIB_USE_LAPACK is #defined then the xGESVD routine
              from LAPACK is used to compute the SVD.
    !*/

// ----------------------------------------------------------------------------------------

    const matrix real_eigenvalues (
        const matrix_exp& m
    );
    /*!
        requires
            - m.nr() == m.nc()
            - matrix_exp::type == float or double
        ensures
            - returns a matrix E such that:
                - E.nr() == m.nr()
                - E.nc() == 1
                - E contains the real part of all eigenvalues of the matrix m.
                  (note that the eigenvalues are not sorted)
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::type det (
        const matrix_exp& m
    );
    /*!
        requires
            - m is a square matrix
        ensures
            - returns the determinant of m
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::type trace (
        const matrix_exp& m
    );
    /*!
        requires
            - m is a square matrix
        ensures
            - returns the trace of m
              (i.e. returns sum(diag(m)))
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::matrix_type chol (
        const matrix_exp& A
    );
    /*!
        requires
            - A is a square matrix
        ensures
            - if (A has a Cholesky Decomposition) then
                - returns the decomposition of A.  That is, returns a matrix L
                  such that L*trans(L) == A.  L will also be lower triangular.
            - else
                - returns a matrix with the same dimensions as A but it 
                  will have a bogus value.  I.e. it won't be a decomposition.
                  In this case the algorithm returns a partial decomposition.
                - You can tell when chol fails by looking at the lower right
                  element of the returned matrix.  If it is 0 then it means
                  A does not have a cholesky decomposition.  

            - If DLIB_USE_LAPACK is defined then the LAPACK routine xPOTRF 
              is used to compute the cholesky decomposition.
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::matrix_type inv_lower_triangular (
        const matrix_exp& A
    );
    /*!
        requires
            - A is a square matrix
        ensures
            - if (A is lower triangular) then
                - returns the inverse of A. 
            - else
                - returns a matrix with the same dimensions as A but it 
                  will have a bogus value.  I.e. it won't be an inverse.
    !*/

// ----------------------------------------------------------------------------------------

    const matrix_exp::matrix_type inv_upper_triangular (
        const matrix_exp& A
    );
    /*!
        requires
            - A is a square matrix
        ensures
            - if (A is upper triangular) then
                - returns the inverse of A. 
            - else
                - returns a matrix with the same dimensions as A but it 
                  will have a bogus value.  I.e. it won't be an inverse.
    !*/

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
//                             Matrix decomposition classes 
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

    template <
        typename matrix_exp_type
        >
    class lu_decomposition
    {
        /*!
            REQUIREMENTS ON matrix_exp_type
                must be some kind of matrix expression as defined in the 
                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
                The matrix type must also contain float or double values.

            WHAT THIS OBJECT REPRESENTS
                This object represents something that can compute an LU 
                decomposition of a real valued matrix.  That is, for any 
                matrix A it computes matrices L, U, and a pivot vector P such 
                that rowm(A,P) == L*U.

                The LU decomposition with pivoting always exists, even if the matrix is
                singular, so the constructor will never fail.  The primary use of the
                LU decomposition is in the solution of square systems of simultaneous
                linear equations.  This will fail if is_singular() returns true (or
                if A is very nearly singular).

                If DLIB_USE_LAPACK is defined then the LAPACK routine xGETRF 
                is used to compute the LU decomposition.
        !*/

    public:

        const static long NR = matrix_exp_type::NR;
        const static long NC = matrix_exp_type::NC;
        typedef typename matrix_exp_type::type type;
        typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
        typedef typename matrix_exp_type::layout_type layout_type;

        typedef matrix<type,0,0,mem_manager_type,layout_type>  matrix_type;
        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
        typedef matrix<long,NR,1,mem_manager_type,layout_type> pivot_column_vector_type;

        template <typename EXP>
        lu_decomposition (
            const matrix_exp<EXP> &A
        );
        /*!
            requires
                - EXP::type == lu_decomposition::type 
                - A.size() > 0
            ensures
                - #nr() == A.nr()
                - #nc() == A.nc()
                - #is_square() == (A.nr() == A.nc())
                - computes the LU factorization of the given A matrix.
        !*/

        bool is_square (
        ) const;
        /*!
            ensures
                - if (the input A matrix was a square matrix) then
                    - returns true
                - else
                    - returns false
        !*/

        bool is_singular (
        ) const;
        /*!
            requires
                - is_square() == true
            ensures
                - if (the input A matrix is singular) then
                    - returns true
                - else
                    - returns false
        !*/

        long nr(
        ) const;
        /*!
            ensures
                - returns the number of rows in the input matrix
        !*/

        long nc(
        ) const;
        /*!
            ensures
                - returns the number of columns in the input matrix
        !*/

        const matrix_type get_l (
        ) const; 
        /*!
            ensures
                - returns the lower triangular L factor of the LU factorization.  
                - L.nr() == nr()
                - L.nc() == min(nr(),nc())
        !*/

        const matrix_type get_u (
        ) const;
        /*!
            ensures
                - returns the upper triangular U factor of the LU factorization.  
                - U.nr() == min(nr(),nc())
                - U.nc() == nc()
        !*/

        const pivot_column_vector_type& get_pivot (
        ) const;
        /*!
            ensures
                - returns the pivot permutation vector.  That is,
                  if A is the input matrix then this function 
                  returns a vector P such that:
                    - rowm(A,P) == get_l()*get_u() 
                    - P.nr() == A.nr()
        !*/

        type det (
        ) const;
        /*!
            requires
                - is_square() == true
            ensures
                - computes and returns the determinant of the input 
                  matrix using LU factors.
        !*/

        template <typename EXP>
        const matrix_type solve (
            const matrix_exp<EXP> &B
        ) const;
        /*!
            requires
                - EXP::type == lu_decomposition::type
                - is_square() == true
                - B.nr() == nr()
            ensures
                - Let A denote the input matrix to this class's constructor.  
                  Then this function solves A*X == B for X and returns X.  
                - Note that if A is singular (or very close to singular) then
                  the X returned by this function won't fit A*X == B very well (if at all).
        !*/

    }; 

// ----------------------------------------------------------------------------------------

    template <
        typename matrix_exp_type
        >
    class cholesky_decomposition
    {
        /*! 
            REQUIREMENTS ON matrix_exp_type
                must be some kind of matrix expression as defined in the 
                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
                The matrix type must also contain float or double values.

            WHAT THIS OBJECT REPRESENTS
                This object represents something that can compute a cholesky 
                decomposition of a real valued matrix.  That is, for any 
                symmetric, positive definite matrix A, it computes a lower 
                triangular matrix L such that A == L*trans(L).
                
                If the matrix is not symmetric or positive definite, the function
                computes only a partial decomposition.  This can be tested with
                the is_spd() flag.
            
                If DLIB_USE_LAPACK is defined then the LAPACK routine xPOTRF 
                is used to compute the cholesky decomposition.
        !*/

    public:

        const static long NR = matrix_exp_type::NR;
        const static long NC = matrix_exp_type::NC;
        typedef typename matrix_exp_type::type type;
        typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
        typedef typename matrix_exp_type::layout_type layout_type;

        typedef typename matrix_exp_type::matrix_type matrix_type;
        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;

        template <typename EXP>
        cholesky_decomposition(
            const matrix_exp<EXP>& A
        );
        /*!
            requires
                - EXP::type == cholesky_decomposition::type 
                - A.size() > 0
                - A.nr() == A.nc() 
                  (i.e. A must be a square matrix)
            ensures
                - if (A is symmetric positive-definite) then
                    - #is_spd() == true 
                    - Constructs a lower triangular matrix L, such that L*trans(L) == A.
                      and #get_l() == L
                - else
                    - #is_spd() == false
        !*/

        bool is_spd(
        ) const;
        /*!
            ensures
                - if (the input matrix was symmetric positive-definite) then
                    - returns true
                - else
                    - returns false
        !*/

        const matrix_type& get_l(
        ) const;
        /*!
            ensures
                - returns the lower triangular factor, L, such that L*trans(L) == A
                  (where A is the input matrix to this class's constructor)
                - Note that if A is not symmetric positive definite or positive semi-definite
                  then the equation L*trans(L) == A won't hold.  
        !*/

        template <typename EXP>
        const matrix solve (
            const matrix_exp<EXP>& B
        ) const;
        /*!
            requires
                - EXP::type == cholesky_decomposition::type
                - B.nr() == get_l().nr()
                  (i.e. the number of rows in B must match the number of rows in the
                  input matrix A)
            ensures
                - Let A denote the input matrix to this class's constructor.  Then 
                  this function solves A*X = B for X and returns X.  
                - Note that if is_spd() == false or A was really close to being
                  non-SPD then the solver will fail to find an accurate solution.
        !*/

    };

// ----------------------------------------------------------------------------------------

    template <
        typename matrix_exp_type
        >
    class qr_decomposition 
    {
        /*! 
            REQUIREMENTS ON matrix_exp_type
                must be some kind of matrix expression as defined in the 
                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
                The matrix type must also contain float or double values.

            WHAT THIS OBJECT REPRESENTS
                This object represents something that can compute a classical
                QR decomposition of an m-by-n real valued matrix A with m >= n.  

                The QR decomposition is an m-by-n orthogonal matrix Q and an 
                n-by-n upper triangular matrix R so that A == Q*R. The QR decomposition 
                always exists, even if the matrix does not have full rank, so the 
                constructor will never fail.  The primary use of the QR decomposition 
                is in the least squares solution of non-square systems of simultaneous 
                linear equations.  This will fail if is_full_rank() returns false or
                A is very nearly not full rank.

                The Q and R factors can be retrieved via the get_q() and get_r()
                methods. Furthermore, a solve() method is provided to find the
                least squares solution of Ax=b using the QR factors.  

                If DLIB_USE_LAPACK is #defined then the xGEQRF routine
                from LAPACK is used to compute the QR decomposition.
        !*/

    public:

        const static long NR = matrix_exp_type::NR;
        const static long NC = matrix_exp_type::NC;
        typedef typename matrix_exp_type::type type;
        typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
        typedef typename matrix_exp_type::layout_type layout_type;

        typedef matrix<type,0,0,mem_manager_type,layout_type> matrix_type;

        template <typename EXP>
        qr_decomposition(
            const matrix_exp<EXP>& A
        );
        /*!
            requires
                - EXP::type == qr_decomposition::type
                - A.nr() >= A.nc()
                - A.size() > 0
            ensures
                - #nr() == A.nr()
                - #nc() == A.nc()
                - computes the QR decomposition of the given A matrix.
        !*/

        bool is_full_rank(
        ) const;
        /*!
            ensures
                - if (the input A matrix had full rank) then
                    - returns true
                - else
                    - returns false
        !*/

        long nr(
        ) const;
        /*!
            ensures
                - returns the number of rows in the input matrix
        !*/

        long nc(
        ) const;
        /*!
            ensures
                - returns the number of columns in the input matrix
        !*/

        const matrix_type get_r (
        ) const;
        /*!
            ensures
                - returns a matrix R such that: 
                    - R is the upper triangular factor, R, of the QR factorization
                    - get_q()*R == input matrix A
                    - R.nr() == nc()
                    - R.nc() == nc()
        !*/

        const matrix_type get_q (
        ) const;
        /*!
            ensures
                - returns a matrix Q such that:
                    - Q is the economy-sized orthogonal factor Q from the QR 
                      factorization.  
                    - trans(Q)*Q == identity matrix
                    - Q*get_r() == input matrix A 
                    - Q.nr() == nr()
                    - Q.nc() == nc()
        !*/

        template <typename EXP>
        const matrix_type solve (
            const matrix_exp<EXP>& B
        ) const;
        /*!
            requires
                - EXP::type == qr_decomposition::type
                - B.nr() == nr()
            ensures
                - Let A denote the input matrix to this class's constructor.  
                  Then this function finds the least squares solution to the equation A*X = B 
                  and returns X.  X has the following properties: 
                    - X is the matrix that minimizes the two norm of A*X-B.  That is, it
                      minimizes sum(squared(A*X - B)).
                    - X.nr() == nc()
                    - X.nc() == B.nc()
                - Note that this function will fail to output a good solution if is_full_rank() == false
                  or the A matrix is close to not being full rank.
        !*/

    };

// ----------------------------------------------------------------------------------------

    template <
        typename matrix_exp_type
        >
    class eigenvalue_decomposition
    {
        /*!
            REQUIREMENTS ON matrix_exp_type
                must be some kind of matrix expression as defined in the 
                dlib/matrix/matrix_abstract.h file.   (e.g. a dlib::matrix object)
                The matrix type must also contain float or double values.

            WHAT THIS OBJECT REPRESENTS
                This object represents something that can compute an eigenvalue 
                decomposition of a real valued matrix.   So it gives 
                you the set of eigenvalues and eigenvectors for a matrix.   

                Let A denote the input matrix to this object's constructor.  Then 
                what this object does is it finds two matrices, D and V, such that
                    - A*V == V*D
                Where V is a square matrix that contains all the eigenvectors
                of the A matrix (each column of V is an eigenvector) and
                D is a diagonal matrix containing the eigenvalues of A.


                It is important to note that if A is symmetric or non-symmetric you
                get somewhat different results.  If A is a symmetric matrix (i.e. A == trans(A))
                then:
                    - All the eigenvalues and eigenvectors of A are real numbers. 
                        - Because of this there isn't really any point in using the
                          part of this class's interface that returns complex matrices.
                          All you need are the get_real_eigenvalues() and
                          get_pseudo_v() functions.  
                    - V*trans(V) should be equal to the identity matrix.  That is, all the
                      eigenvectors in V should be orthonormal. 
                        - So A == V*D*trans(V)
                    - If DLIB_USE_LAPACK is #defined then this object uses the xSYEVR LAPACK
                      routine.

                On the other hand, if A is not symmetric then:
                    - Some of the eigenvalues and eigenvectors might be complex numbers.  
                        - An eigenvalue is complex if and only if its corresponding eigenvector 
                          is complex.  So you can check for this case by just checking 
                          get_imag_eigenvalues() to see if any values are non-zero.  You don't 
                          have to check the V matrix as well.
                    - V*trans(V) won't be equal to the identity matrix but it is usually
                      invertible.  So A == V*D*inv(V) is usually a valid statement but
                      A == V*D*trans(V) won't be.
                    - If DLIB_USE_LAPACK is #defined then this object uses the xGEEV LAPACK
                      routine.
        !*/

    public:

        const static long NR = matrix_exp_type::NR;
        const static long NC = matrix_exp_type::NC;
        typedef typename matrix_exp_type::type type;
        typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
        typedef typename matrix_exp_type::layout_type layout_type;

        typedef typename matrix_exp_type::matrix_type matrix_type;
        typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;

        typedef matrix<std::complex<type>,0,0,mem_manager_type,layout_type> complex_matrix_type;
        typedef matrix<std::complex<type>,NR,1,mem_manager_type,layout_type> complex_column_vector_type;


        template <typename EXP>
        eigenvalue_decomposition(
            const matrix_exp<EXP>& A
        ); 
        /*!
            requires
                - A.nr() == A.nc() 
                - A.size() > 0
                - EXP::type == eigenvalue_decomposition::type 
            ensures
                - #dim() == A.nr()
                - computes the eigenvalue decomposition of A.  
                - #get_eigenvalues() == the eigenvalues of A
                - #get_v() == all the eigenvectors of A
        !*/

        template <typename EXP>
        eigenvalue_decomposition(
            const matrix_op<op_make_symmetric<EXP> >& A
        ); 
        /*!
            requires
                - A.nr() == A.nc() 
                - A.size() > 0
                - EXP::type == eigenvalue_decomposition::type 
            ensures
                - #dim() == A.nr()
                - computes the eigenvalue decomposition of the symmetric matrix A.  Does so
                  using a method optimized for symmetric matrices.
                - #get_eigenvalues() == the eigenvalues of A
                - #get_v() == all the eigenvectors of A
                - moreover, since A is symmetric there won't be any imaginary eigenvalues. So 
                  we will have:
                    - #get_imag_eigenvalues() == 0
                    - #get_real_eigenvalues() == the eigenvalues of A
                    - #get_pseudo_v() == all the eigenvectors of A
                    - diagm(#get_real_eigenvalues()) == #get_pseudo_d()

                Note that the symmetric matrix operator is created by the
                dlib::make_symmetric() function.  This function simply reflects
                the lower triangular part of a square matrix into the upper triangular
                part to create a symmetric matrix.  It can also be used to denote that a 
                matrix is already symmetric using the C++ type system.
        !*/

        long dim (
        ) const;
        /*!
            ensures
                - dim() == the number of rows/columns in the input matrix A 
        !*/

        const complex_column_vector_type get_eigenvalues (
        ) const;
        /*!
            ensures
                - returns diag(get_d()).  That is, returns a 
                  vector that contains the eigenvalues of the input 
                  matrix.
                - the returned vector has dim() rows
                - the eigenvalues are not sorted in any particular way
        !*/

        const column_vector_type& get_real_eigenvalues (
        ) const;
        /*! 
            ensures
                - returns the real parts of the eigenvalues.  That is,
                  returns real(get_eigenvalues()) 
                - the returned vector has dim() rows
                - the eigenvalues are not sorted in any particular way
        !*/

        const column_vector_type& get_imag_eigenvalues (
        ) const;
        /*! 
            ensures
                - returns the imaginary parts of the eigenvalues.  That is,
                  returns imag(get_eigenvalues()) 
                - the returned vector has dim() rows
                - the eigenvalues are not sorted in any particular way
        !*/

        const complex_matrix_type get_v (
        ) const;
        /*!
            ensures
                - returns the eigenvector matrix V that is 
                  dim() rows by dim() columns
                - Each column in V is one of the eigenvectors of the input 
                  matrix
        !*/

        const complex_matrix_type get_d (
        ) const; 
        /*!
            ensures
                - returns a matrix D such that:
                    - D.nr() == dim()
                    - D.nc() == dim()
                    - diag(D) == get_eigenvalues()
                      (i.e. the diagonal of D contains all the eigenvalues in the input matrix)
                    - all off diagonal elements of D are set to 0
        !*/

        const matrix_type& get_pseudo_v (
        ) const;
        /*!
            ensures
                - returns a matrix that is dim() rows by dim() columns
                - Let A denote the input matrix given to this object's constructor.
                - if (A has any imaginary eigenvalues) then
                    - returns the pseudo-eigenvector matrix V  
                    - The matrix V returned by this function is structured such that:
                        - A*V == V*get_pseudo_d()
                - else
                    - returns the eigenvector matrix V with A's eigenvectors as
                      the columns of V
                    - A*V == V*diagm(get_real_eigenvalues())
        !*/

        const matrix_type get_pseudo_d (
        ) const; 
        /*!
            ensures
                - The returned matrix is dim() rows by dim() columns
                - Computes and returns the block diagonal eigenvalue matrix.
                  If the original matrix A is not symmetric, then the eigenvalue 
                  matrix D is block diagonal with the real eigenvalues in 1-by-1 
                  blocks and any complex eigenvalues,
                  a + i*b, in 2-by-2 blocks, (a, b; -b, a).  That is, if the complex
                  eigenvalues look like

                      u + iv     .        .          .      .    .
                        .      u - iv     .          .      .    .
                        .        .      a + ib       .      .    .
                        .        .        .        a - ib   .    .
                        .        .        .          .      x    .
                        .        .        .          .      .    y

                  Then D looks like

                        u        v        .          .      .    .
                       -v        u        .          .      .    . 
                        .        .        a          b      .    .
                        .        .       -b          a      .    .
                        .        .        .          .      x    .
                        .        .        .          .      .    y

                  This keeps V (The V you get from get_pseudo_v()) a real matrix in both 
                  symmetric and non-symmetric cases, and A*V = V*D.
                - the eigenvalues are not sorted in any particular way
        !*/

    };

// ----------------------------------------------------------------------------------------

}

#endif // DLIB_MATRIx_LA_FUNCTS_ABSTRACT_