//  $Id: mmdb_math_linalg.h $
//  =================================================================
//
//   CCP4 Coordinate Library: support of coordinate-related
//   functionality in protein crystallography applications.
//
//   Copyright (C) Eugene Krissinel 2000-2013.
//
//    This library is free software: you can redistribute it and/or
//    modify it under the terms of the GNU Lesser General Public
//    License version 3, modified in accordance with the provisions
//    of the license to address the requirements of UK law.
//
//    You should have received a copy of the modified GNU Lesser
//    General Public License along with this library. If not, copies
//    may be downloaded from http://www.ccp4.ac.uk/ccp4license.php
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU Lesser General Public License for more details.
//
//  =================================================================
//
//    12.09.13   <--  Date of Last Modification.
//                   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//  -----------------------------------------------------------------
//
//  **** Module  :  LinAlg  <interface>
//       ~~~~~~~~~
//  **** Project :  MMDB  ( MacroMolecular Data Base )
//       ~~~~~~~~~
//
//  (C) E.Krissinel  2000-2013
//
//  =================================================================
//
//

#ifndef __MMDB_MATH_LinAlg__
#define __MMDB_MATH_LinAlg__

#include "mmdb_mattype.h"

namespace mmdb  {

  namespace math  {

    //  ==========================  Jacobi  =============================


    ///  Diagonalization of symmetric matrices A[1..N][1..N]
    /// by the method of Jacobi.
    extern void  Jacobi ( int     N,     //!< dimension of the matrix
                          rmatrix A,     //!< matrix to diagonalize; the
                                         /// lower triangle, except the
                                         /// diagonal, will remain unchanged
                          rmatrix T,     //!< eigenvectors placed as columns
                          rvector Eigen, //!< vector of eigenvalues, orderd
                                         /// by increasing
                          rvector Aik,   //!< working array
                          int &  Signal  //!< 0 <=> Ok, ItMax <=> iteration
                                         /// limit exchausted.
                        );


    //  A5.5.2  :  Perturbed Cholesky Decomposition
    extern void  PbCholDecomp ( int        N,
                                rvector    HDiag,
                                realtype   MaxOff,
                                realtype   MachEps,
                                rmatrix    L,
                                realtype & MaxAdd );

    //  A3.2.3a  :  Cholesky's   L - Solution  of
    //              L*Y  =  B  ( given  B )
    extern void  LSolve ( int N, rmatrix L, rvector B, rvector Y );

    //  A3.2.3b  :  Cholesky's   LT - Solution  of
    //              LT*X  =  Y  ( given  Y )
    extern void  LTSolve ( int N, rmatrix L, rvector Y, rvector X );

    //  A3.2.3   :  Solution of the equation    L*LT*S = G
    //              by the  Cholesky's  method
    extern void  ChSolve ( int N, rmatrix L, rvector G, rvector S );


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

    extern void  FastInverse (  int N, rmatrix A, ivector J0,
    //#D                          realtype &  Det,
                                int & Signal );
    //
    //      13.09.90  <--  Last Modification Date
    //                    ------------------------
    //
    // ================================================
    //
    //        Fast Inversion of the matrix  A
    //      by the method of  GAUSS - JORDAN  .
    //
    // ------------------------------------------------
    //
    //          Input  parameters  are  :
    //
    //     N   -   dimension of the matrix
    //     A   -   the matrix [1..N][1..N] to be inverted.
    // ------------------------------------------------
    //
    //     J0  -   integer vector [1..N] for temporal storage
    //
    // ------------------------------------------------
    //
    //          Output parameters  are  :
    //
    //     A   -   the inverted matrix
    //     Signal - the error key :
    //            = 0   <=>   O'K
    //             else
    //            degeneration was found, and
    //            the rang of matrix is  Signal-1.
    //
    //        Variable  Det  may return the determinant
    //     of matrix A.  To obtain it, remove all comments
    //     of form  //#D.
    //
    // ================================================


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

    extern void  SVD ( int    NA,  int     M,    int N,
                       rmatrix A,  rmatrix U,    rmatrix V,
                       rvector W,  rvector RV1,
                       bool MatU,  bool   MatV,
                       int & RetCode );
    //
    //      13.12.01  <--  Last Modification Date
    //                    ------------------------
    //
    // ================================================
    //
    //         The    Singular Value Decomposition
    //    of the matrix  A  by the algorithm from
    //      G.Forsait, M.Malkolm, K.Mouler.  Numerical
    //    methods of mathematical calculations //
    //    M., Mir, 1980.
    //
    //         Matrix  A  is represented as
    //
    //         A  =  U * W * VT
    //
    // ------------------------------------------------
    //
    //  All dimensions are indexed from 1 on.
    //
    // ------------------------------------------------
    //
    //         Input  parameters:
    //
    //     NA  -   number of lines in A. NA may be
    //           equal to M or N  only.  If NA=M
    //           then usual SVD will be made. If MA=N
    //           then matrix A is transposed before
    //           the decomposition, and the meaning of
    //           output parameters  U  and  V  is
    //           swapped (U accepts VT and VT accepts U).
    //           In other words, matrix  A  has physical
    //           dimension of  M x N , same as U and V;
    //           however the logical dimension of it
    //           remains that of  N x M .
    //     M   -   number of lines in  U
    //     N   -   number of columns in  U,V and length
    //           of  W,RV1 .  Always provide  M >= N  !
    //     A   -   matrix [1..M][1..N] or [1..N][1..M]
    //           to be decomposed. The matrix does not
    //           change,  and it may coincide with  U  or
    //           V, if NA=M (in which case  A  does change)
    //     MatU -  compute  U , if set True
    //     MatV -  compute  V , if set True
    //     RV1  -  temporary array [1..N].
    //     U    - should be always supplied as an array of
    //            [1..M][1..N], M>=N .
    //     V    - should be suuplied as an array of
    //            [1..N][1..N] if MatV is True .
    //
    // ------------------------------------------------
    //
    //          Output parameters  are  :
    //
    //     W   -   N non-ordered singular values,
    //           if  RetCode=0. If RetCode<>0, the
    //           RetCode+1 ... N -th values are still
    //           valid
    //     U   -   matrix of right singular vectors
    //           (arranged in columns),  corresponding
    //           to the singular values in  W,  if
    //           RetCode=0 and MatU is True.  If MatU
    //           is False, U  is still used as a
    //           temporary array. If RetCode<>0 then
    //           the  RetCode+1 ... N -th vectors
    //           are  valid
    //     V   -   matrix of left singular vectors
    //           (arranged in columns),  corresponding
    //           to the singular values in  W,  if
    //           RetCode=0 and MatV is True. If MatV
    //           is False, V is not used and may be set
    //           to NULL. If RetCode<>0 then the
    //           RetCode+1 ... N -th vectors are valid
    //     RetCode - the error key :
    //            = 0   <=>   O'K
    //              else
    //            = k, if the k-th singular value
    //                 was not computed after 30 iterations.
    //
    // ------------------------------------------------
    //
    //          Key  Variables  are  :
    //
    //     ItnLimit  -  the limit for iterations
    //
    //     This routine does not use any machine-dependent
    //  constants.
    //
    // ================================================
    //
    //

    extern void  OrderSVD ( int M, int N, rmatrix U, rmatrix V,
                            rvector W, bool MatU, bool MatV );


  }
}

#endif