| ------------------------------------------------------------------------------ |
| -- -- |
| -- GNAT RUN-TIME COMPONENTS -- |
| -- -- |
| -- SYSTEM.GENERIC_REAL_LAPACK -- |
| -- -- |
| -- B o d y -- |
| -- -- |
| -- Copyright (C) 2006, Free Software Foundation, Inc. -- |
| -- -- |
| -- GNAT is free software; you can redistribute it and/or modify it under -- |
| -- terms of the GNU General Public License as published by the Free Soft- -- |
| -- ware Foundation; either version 2, or (at your option) any later ver- -- |
| -- sion. GNAT is distributed in the hope that it will be useful, but WITH- -- |
| -- OUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- |
| -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -- |
| -- for more details. You should have received a copy of the GNU General -- |
| -- Public License distributed with GNAT; see file COPYING. If not, write -- |
| -- to the Free Software Foundation, 51 Franklin Street, Fifth Floor, -- |
| -- Boston, MA 02110-1301, USA. -- |
| -- -- |
| -- As a special exception, if other files instantiate generics from this -- |
| -- unit, or you link this unit with other files to produce an executable, -- |
| -- this unit does not by itself cause the resulting executable to be -- |
| -- covered by the GNU General Public License. This exception does not -- |
| -- however invalidate any other reasons why the executable file might be -- |
| -- covered by the GNU Public License. -- |
| -- -- |
| -- GNAT was originally developed by the GNAT team at New York University. -- |
| -- Extensive contributions were provided by Ada Core Technologies Inc. -- |
| -- -- |
| ------------------------------------------------------------------------------ |
| |
| with Ada.Unchecked_Conversion; use Ada; |
| with Interfaces; use Interfaces; |
| with Interfaces.Fortran; use Interfaces.Fortran; |
| with Interfaces.Fortran.BLAS; use Interfaces.Fortran.BLAS; |
| with Interfaces.Fortran.LAPACK; use Interfaces.Fortran.LAPACK; |
| with System.Generic_Array_Operations; use System.Generic_Array_Operations; |
| |
| package body System.Generic_Real_LAPACK is |
| |
| Is_Real : constant Boolean := |
| Real'Machine_Mantissa = Fortran.Real'Machine_Mantissa |
| and then Fortran.Real (Real'First) = Fortran.Real'First |
| and then Fortran.Real (Real'Last) = Fortran.Real'Last; |
| |
| Is_Double_Precision : constant Boolean := |
| Real'Machine_Mantissa = |
| Double_Precision'Machine_Mantissa |
| and then |
| Double_Precision (Real'First) = |
| Double_Precision'First |
| and then |
| Double_Precision (Real'Last) = |
| Double_Precision'Last; |
| |
| -- Local subprograms |
| |
| function To_Double_Precision (X : Real) return Double_Precision; |
| pragma Inline_Always (To_Double_Precision); |
| |
| function To_Real (X : Double_Precision) return Real; |
| pragma Inline_Always (To_Real); |
| |
| -- Instantiations |
| |
| function To_Double_Precision is new |
| Vector_Elementwise_Operation |
| (X_Scalar => Real, |
| Result_Scalar => Double_Precision, |
| X_Vector => Real_Vector, |
| Result_Vector => Double_Precision_Vector, |
| Operation => To_Double_Precision); |
| |
| function To_Real is new |
| Vector_Elementwise_Operation |
| (X_Scalar => Double_Precision, |
| Result_Scalar => Real, |
| X_Vector => Double_Precision_Vector, |
| Result_Vector => Real_Vector, |
| Operation => To_Real); |
| |
| function To_Double_Precision is new |
| Matrix_Elementwise_Operation |
| (X_Scalar => Real, |
| Result_Scalar => Double_Precision, |
| X_Matrix => Real_Matrix, |
| Result_Matrix => Double_Precision_Matrix, |
| Operation => To_Double_Precision); |
| |
| function To_Real is new |
| Matrix_Elementwise_Operation |
| (X_Scalar => Double_Precision, |
| Result_Scalar => Real, |
| X_Matrix => Double_Precision_Matrix, |
| Result_Matrix => Real_Matrix, |
| Operation => To_Real); |
| |
| function To_Double_Precision (X : Real) return Double_Precision is |
| begin |
| return Double_Precision (X); |
| end To_Double_Precision; |
| |
| function To_Real (X : Double_Precision) return Real is |
| begin |
| return Real (X); |
| end To_Real; |
| |
| ----------- |
| -- getrf -- |
| ----------- |
| |
| procedure getrf |
| (M : Natural; |
| N : Natural; |
| A : in out Real_Matrix; |
| Ld_A : Positive; |
| I_Piv : out Integer_Vector; |
| Info : access Integer) |
| is |
| begin |
| if Is_Real then |
| declare |
| type A_Ptr is |
| access all BLAS.Real_Matrix (A'Range (1), A'Range (2)); |
| function Conv_A is new Unchecked_Conversion (Address, A_Ptr); |
| begin |
| sgetrf (M, N, Conv_A (A'Address).all, Ld_A, |
| LAPACK.Integer_Vector (I_Piv), Info); |
| end; |
| |
| elsif Is_Double_Precision then |
| declare |
| type A_Ptr is |
| access all Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| function Conv_A is new Unchecked_Conversion (Address, A_Ptr); |
| begin |
| dgetrf (M, N, Conv_A (A'Address).all, Ld_A, |
| LAPACK.Integer_Vector (I_Piv), Info); |
| end; |
| |
| else |
| declare |
| DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| begin |
| DP_A := To_Double_Precision (A); |
| dgetrf (M, N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), Info); |
| A := To_Real (DP_A); |
| end; |
| end if; |
| end getrf; |
| |
| ----------- |
| -- getri -- |
| ----------- |
| |
| procedure getri |
| (N : Natural; |
| A : in out Real_Matrix; |
| Ld_A : Positive; |
| I_Piv : Integer_Vector; |
| Work : in out Real_Vector; |
| L_Work : Integer; |
| Info : access Integer) |
| is |
| begin |
| if Is_Real then |
| declare |
| type A_Ptr is |
| access all BLAS.Real_Matrix (A'Range (1), A'Range (2)); |
| type Work_Ptr is |
| access all BLAS.Real_Vector (Work'Range); |
| function Conv_A is new Unchecked_Conversion (Address, A_Ptr); |
| function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); |
| begin |
| sgetri (N, Conv_A (A'Address).all, Ld_A, |
| LAPACK.Integer_Vector (I_Piv), |
| Conv_Work (Work'Address).all, L_Work, |
| Info); |
| end; |
| |
| elsif Is_Double_Precision then |
| declare |
| type A_Ptr is |
| access all Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| type Work_Ptr is |
| access all Double_Precision_Vector (Work'Range); |
| function Conv_A is new Unchecked_Conversion (Address, A_Ptr); |
| function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); |
| begin |
| dgetri (N, Conv_A (A'Address).all, Ld_A, |
| LAPACK.Integer_Vector (I_Piv), |
| Conv_Work (Work'Address).all, L_Work, |
| Info); |
| end; |
| |
| else |
| declare |
| DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| DP_Work : Double_Precision_Vector (Work'Range); |
| begin |
| DP_A := To_Double_Precision (A); |
| dgetri (N, DP_A, Ld_A, LAPACK.Integer_Vector (I_Piv), |
| DP_Work, L_Work, Info); |
| A := To_Real (DP_A); |
| Work (1) := To_Real (DP_Work (1)); |
| end; |
| end if; |
| end getri; |
| |
| ----------- |
| -- getrs -- |
| ----------- |
| |
| procedure getrs |
| (Trans : access constant Character; |
| N : Natural; |
| N_Rhs : Natural; |
| A : Real_Matrix; |
| Ld_A : Positive; |
| I_Piv : Integer_Vector; |
| B : in out Real_Matrix; |
| Ld_B : Positive; |
| Info : access Integer) |
| is |
| begin |
| if Is_Real then |
| declare |
| subtype A_Type is BLAS.Real_Matrix (A'Range (1), A'Range (2)); |
| type B_Ptr is |
| access all BLAS.Real_Matrix (B'Range (1), B'Range (2)); |
| function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type); |
| function Conv_B is new Unchecked_Conversion (Address, B_Ptr); |
| begin |
| sgetrs (Trans, N, N_Rhs, |
| Conv_A (A), Ld_A, |
| LAPACK.Integer_Vector (I_Piv), |
| Conv_B (B'Address).all, Ld_B, |
| Info); |
| end; |
| |
| elsif Is_Double_Precision then |
| declare |
| subtype A_Type is |
| Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| type B_Ptr is |
| access all Double_Precision_Matrix (B'Range (1), B'Range (2)); |
| function Conv_A is new Unchecked_Conversion (Real_Matrix, A_Type); |
| function Conv_B is new Unchecked_Conversion (Address, B_Ptr); |
| begin |
| dgetrs (Trans, N, N_Rhs, |
| Conv_A (A), Ld_A, |
| LAPACK.Integer_Vector (I_Piv), |
| Conv_B (B'Address).all, Ld_B, |
| Info); |
| end; |
| |
| else |
| declare |
| DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| DP_B : Double_Precision_Matrix (B'Range (1), B'Range (2)); |
| begin |
| DP_A := To_Double_Precision (A); |
| DP_B := To_Double_Precision (B); |
| dgetrs (Trans, N, N_Rhs, |
| DP_A, Ld_A, |
| LAPACK.Integer_Vector (I_Piv), |
| DP_B, Ld_B, |
| Info); |
| B := To_Real (DP_B); |
| end; |
| end if; |
| end getrs; |
| |
| ----------- |
| -- orgtr -- |
| ----------- |
| |
| procedure orgtr |
| (Uplo : access constant Character; |
| N : Natural; |
| A : in out Real_Matrix; |
| Ld_A : Positive; |
| Tau : Real_Vector; |
| Work : out Real_Vector; |
| L_Work : Integer; |
| Info : access Integer) |
| is |
| begin |
| if Is_Real then |
| declare |
| type A_Ptr is |
| access all BLAS.Real_Matrix (A'Range (1), A'Range (2)); |
| subtype Tau_Type is BLAS.Real_Vector (Tau'Range); |
| type Work_Ptr is |
| access all BLAS.Real_Vector (Work'Range); |
| function Conv_A is new Unchecked_Conversion (Address, A_Ptr); |
| function Conv_Tau is |
| new Unchecked_Conversion (Real_Vector, Tau_Type); |
| function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); |
| begin |
| sorgtr (Uplo, N, |
| Conv_A (A'Address).all, Ld_A, |
| Conv_Tau (Tau), |
| Conv_Work (Work'Address).all, L_Work, |
| Info); |
| end; |
| |
| elsif Is_Double_Precision then |
| declare |
| type A_Ptr is |
| access all Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| subtype Tau_Type is Double_Precision_Vector (Tau'Range); |
| type Work_Ptr is |
| access all Double_Precision_Vector (Work'Range); |
| function Conv_A is new Unchecked_Conversion (Address, A_Ptr); |
| function Conv_Tau is |
| new Unchecked_Conversion (Real_Vector, Tau_Type); |
| function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); |
| begin |
| dorgtr (Uplo, N, |
| Conv_A (A'Address).all, Ld_A, |
| Conv_Tau (Tau), |
| Conv_Work (Work'Address).all, L_Work, |
| Info); |
| end; |
| |
| else |
| declare |
| DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| DP_Work : Double_Precision_Vector (Work'Range); |
| DP_Tau : Double_Precision_Vector (Tau'Range); |
| begin |
| DP_A := To_Double_Precision (A); |
| DP_Tau := To_Double_Precision (Tau); |
| dorgtr (Uplo, N, DP_A, Ld_A, DP_Tau, DP_Work, L_Work, Info); |
| A := To_Real (DP_A); |
| Work (1) := To_Real (DP_Work (1)); |
| end; |
| end if; |
| end orgtr; |
| |
| ----------- |
| -- steqr -- |
| ----------- |
| |
| procedure steqr |
| (Comp_Z : access constant Character; |
| N : Natural; |
| D : in out Real_Vector; |
| E : in out Real_Vector; |
| Z : in out Real_Matrix; |
| Ld_Z : Positive; |
| Work : out Real_Vector; |
| Info : access Integer) |
| is |
| begin |
| if Is_Real then |
| declare |
| type D_Ptr is access all BLAS.Real_Vector (D'Range); |
| type E_Ptr is access all BLAS.Real_Vector (E'Range); |
| type Z_Ptr is |
| access all BLAS.Real_Matrix (Z'Range (1), Z'Range (2)); |
| type Work_Ptr is |
| access all BLAS.Real_Vector (Work'Range); |
| function Conv_D is new Unchecked_Conversion (Address, D_Ptr); |
| function Conv_E is new Unchecked_Conversion (Address, E_Ptr); |
| function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr); |
| function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); |
| begin |
| ssteqr (Comp_Z, N, |
| Conv_D (D'Address).all, |
| Conv_E (E'Address).all, |
| Conv_Z (Z'Address).all, |
| Ld_Z, |
| Conv_Work (Work'Address).all, |
| Info); |
| end; |
| |
| elsif Is_Double_Precision then |
| declare |
| type D_Ptr is access all Double_Precision_Vector (D'Range); |
| type E_Ptr is access all Double_Precision_Vector (E'Range); |
| type Z_Ptr is |
| access all Double_Precision_Matrix (Z'Range (1), Z'Range (2)); |
| type Work_Ptr is |
| access all Double_Precision_Vector (Work'Range); |
| function Conv_D is new Unchecked_Conversion (Address, D_Ptr); |
| function Conv_E is new Unchecked_Conversion (Address, E_Ptr); |
| function Conv_Z is new Unchecked_Conversion (Address, Z_Ptr); |
| function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); |
| begin |
| dsteqr (Comp_Z, N, |
| Conv_D (D'Address).all, |
| Conv_E (E'Address).all, |
| Conv_Z (Z'Address).all, |
| Ld_Z, |
| Conv_Work (Work'Address).all, |
| Info); |
| end; |
| |
| else |
| declare |
| DP_D : Double_Precision_Vector (D'Range); |
| DP_E : Double_Precision_Vector (E'Range); |
| DP_Z : Double_Precision_Matrix (Z'Range (1), Z'Range (2)); |
| DP_Work : Double_Precision_Vector (Work'Range); |
| begin |
| DP_D := To_Double_Precision (D); |
| DP_E := To_Double_Precision (E); |
| |
| if Comp_Z.all = 'V' then |
| DP_Z := To_Double_Precision (Z); |
| end if; |
| |
| dsteqr (Comp_Z, N, DP_D, DP_E, DP_Z, Ld_Z, DP_Work, Info); |
| |
| D := To_Real (DP_D); |
| E := To_Real (DP_E); |
| Z := To_Real (DP_Z); |
| end; |
| end if; |
| end steqr; |
| |
| ----------- |
| -- sterf -- |
| ----------- |
| |
| procedure sterf |
| (N : Natural; |
| D : in out Real_Vector; |
| E : in out Real_Vector; |
| Info : access Integer) |
| is |
| begin |
| if Is_Real then |
| declare |
| type D_Ptr is access all BLAS.Real_Vector (D'Range); |
| type E_Ptr is access all BLAS.Real_Vector (E'Range); |
| function Conv_D is new Unchecked_Conversion (Address, D_Ptr); |
| function Conv_E is new Unchecked_Conversion (Address, E_Ptr); |
| begin |
| ssterf (N, Conv_D (D'Address).all, Conv_E (E'Address).all, Info); |
| end; |
| |
| elsif Is_Double_Precision then |
| declare |
| type D_Ptr is access all Double_Precision_Vector (D'Range); |
| type E_Ptr is access all Double_Precision_Vector (E'Range); |
| function Conv_D is new Unchecked_Conversion (Address, D_Ptr); |
| function Conv_E is new Unchecked_Conversion (Address, E_Ptr); |
| begin |
| dsterf (N, Conv_D (D'Address).all, Conv_E (E'Address).all, Info); |
| end; |
| |
| else |
| declare |
| DP_D : Double_Precision_Vector (D'Range); |
| DP_E : Double_Precision_Vector (E'Range); |
| |
| begin |
| DP_D := To_Double_Precision (D); |
| DP_E := To_Double_Precision (E); |
| |
| dsterf (N, DP_D, DP_E, Info); |
| |
| D := To_Real (DP_D); |
| E := To_Real (DP_E); |
| end; |
| end if; |
| end sterf; |
| |
| ----------- |
| -- sytrd -- |
| ----------- |
| |
| procedure sytrd |
| (Uplo : access constant Character; |
| N : Natural; |
| A : in out Real_Matrix; |
| Ld_A : Positive; |
| D : out Real_Vector; |
| E : out Real_Vector; |
| Tau : out Real_Vector; |
| Work : out Real_Vector; |
| L_Work : Integer; |
| Info : access Integer) |
| is |
| begin |
| if Is_Real then |
| declare |
| type A_Ptr is |
| access all BLAS.Real_Matrix (A'Range (1), A'Range (2)); |
| type D_Ptr is access all BLAS.Real_Vector (D'Range); |
| type E_Ptr is access all BLAS.Real_Vector (E'Range); |
| type Tau_Ptr is access all BLAS.Real_Vector (Tau'Range); |
| type Work_Ptr is |
| access all BLAS.Real_Vector (Work'Range); |
| function Conv_A is new Unchecked_Conversion (Address, A_Ptr); |
| function Conv_D is new Unchecked_Conversion (Address, D_Ptr); |
| function Conv_E is new Unchecked_Conversion (Address, E_Ptr); |
| function Conv_Tau is new Unchecked_Conversion (Address, Tau_Ptr); |
| function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); |
| begin |
| ssytrd (Uplo, N, |
| Conv_A (A'Address).all, Ld_A, |
| Conv_D (D'Address).all, |
| Conv_E (E'Address).all, |
| Conv_Tau (Tau'Address).all, |
| Conv_Work (Work'Address).all, |
| L_Work, |
| Info); |
| end; |
| |
| elsif Is_Double_Precision then |
| declare |
| type A_Ptr is |
| access all Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| type D_Ptr is access all Double_Precision_Vector (D'Range); |
| type E_Ptr is access all Double_Precision_Vector (E'Range); |
| type Tau_Ptr is access all Double_Precision_Vector (Tau'Range); |
| type Work_Ptr is |
| access all Double_Precision_Vector (Work'Range); |
| function Conv_A is new Unchecked_Conversion (Address, A_Ptr); |
| function Conv_D is new Unchecked_Conversion (Address, D_Ptr); |
| function Conv_E is new Unchecked_Conversion (Address, E_Ptr); |
| function Conv_Tau is new Unchecked_Conversion (Address, Tau_Ptr); |
| function Conv_Work is new Unchecked_Conversion (Address, Work_Ptr); |
| begin |
| dsytrd (Uplo, N, |
| Conv_A (A'Address).all, Ld_A, |
| Conv_D (D'Address).all, |
| Conv_E (E'Address).all, |
| Conv_Tau (Tau'Address).all, |
| Conv_Work (Work'Address).all, |
| L_Work, |
| Info); |
| end; |
| |
| else |
| declare |
| DP_A : Double_Precision_Matrix (A'Range (1), A'Range (2)); |
| DP_D : Double_Precision_Vector (D'Range); |
| DP_E : Double_Precision_Vector (E'Range); |
| DP_Tau : Double_Precision_Vector (Tau'Range); |
| DP_Work : Double_Precision_Vector (Work'Range); |
| begin |
| DP_A := To_Double_Precision (A); |
| |
| dsytrd (Uplo, N, DP_A, Ld_A, DP_D, DP_E, DP_Tau, |
| DP_Work, L_Work, Info); |
| |
| if L_Work /= -1 then |
| A := To_Real (DP_A); |
| D := To_Real (DP_D); |
| E := To_Real (DP_E); |
| Tau := To_Real (DP_Tau); |
| end if; |
| |
| Work (1) := To_Real (DP_Work (1)); |
| end; |
| end if; |
| end sytrd; |
| |
| end System.Generic_Real_LAPACK; |