2.46 Generate the adjugate matrix for square matrix

Given square matrix such as

\[\begin {pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ \end {pmatrix} \]

find the adjugate matrix which is

\[\begin {pmatrix} -3 & 6 & -3 \\ 6 & -12 & 6 \\ -3 & 6 & -3 \\ \end {pmatrix} \]

Mathematica <<"Combinatorica`" (*not needed in V9*) mat = {{1,2,3},{4,5,6},{7,8,9}}; cof = Table[Cofactor[mat,{i,j}], {i,3},{j,3}]; adjugate = Transpose[cof]

{ {-3, 6, -3}, { 6, -12, 6}, {-3, 6, -3}}

 

Maple A:=Matrix([[1,2,3],[4,5,6],[7,8,9]]); LinearAlgebra:-Adjoint(A)

\[ \left [\begin {array}{ccc} -3 & 6 & -3 \\ 6 & -12 & 6 \\ -3 & 6 & -3 \end {array}\right ] \]

 

Matlab Will try to find function in Matlab to do this. But for non-singular matrices only the direct method of finding the inverse and then multiplying by the determinant to recover the adjunct can be used.

A=[1 2 3;4 5 6;7 8 9] det(A)*inv(A)	Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.541976e-18. ans = -3.0000 6.0000 -3.0000 6.0000 -12.0000 6.0000 -3.0000 6.0000 -3.0000
The following is due to Matt J from the matlab newsgroup A=[1 2 3; 4 5 6; 7 8 9] c=poly(A); adjunct_A=polyvalm(c(1:end-1),A)	adjunct_A = -3.0000 6.0000 -3.0000 6.0000 -12.0000 6.0000 -3.0000 6.0000 -3.0000

 

Fortran

Thanks goes to James Van Buskirk and Louisa from comp.lang.fortran for the review and suggestions which lead to improvements of the code.

!-- Find the Adjugate of a square matrix

!-- Version 6/22/2012 by Nasser M. Abbasi gfortran 4.6.3 
program t44 
implicit none 
integer, parameter :: n=3 
integer :: INFO,i,j,ii,IPIV(n-1),number_row_exchange 
real (kind=kind(0.0d0)) :: A(n,n),B(n-1,n-1),C(n,n) 
logical :: M(n,n) 
 
A(1,:) = [1, 2, 3]; 
A(2,:) = [4, 5, 6]; 
A(3,:) = [7, 8, 9]; 
 
DO j=1,n 
   DO i=1,n 
 
      M = .true. 
      M(:,j) = .false. 
      M(i,:) = .false. 
      B = reshape(pack(A, M),[n-1,n-1]) 
 
      !-- LU decomposition in order to obtain determinant 
      CALL DGETRF( n-1, n-1, B, n-1, IPIV, INFO ) 
 
      !-- determinant of U is the product of diagonal 
      C(i,j)= PRODUCT( [(B(ii,ii),ii=1,n-1)] ) 
 
      !-- Adjust sign based on number of row exchanges and (-1)^(i+j) 
      number_row_exchange = count(IPIV /= [(ii,ii=1,n-1)]) 
      C(i,j) = C(i,j)*(1-2*MODULO(number_row_exchange+i+j,2)) 
 
   END DO 
END DO 
write(*,'(3F15.4)') C 
end program t44

#compile, link and run 
export LD_LIBRARY_PATH=/usr/lib/atlas-base/:/usr/lib/atlas-base/atlas/ 
gfortran -fcheck=all -Wall t44.f90 -L/usr/lib/atlas-base/atlas/ -lblas -llapack 
>./a.out 
        -3.0000         6.0000        -3.0000 
         6.0000       -12.0000         6.0000 
        -3.0000         6.0000        -3.0000 
!>dpkg -l|grep -E "(blas|lapack)" 
!ii  libblas3gf        1.2.20110419-2ubuntu1  Basic Linear Algebra Reference implementations, shared library 
!ii  liblapack-dev     3.3.1-1                library of linear algebra routines 3 - static version 
!ii  liblapack-doc     3.3.1-1                library of linear algebra routines 3 - documentation 
!ii  liblapack3gf      3.3.1-1                library of linear algebra routines 3 - shared version 
!ii  libopenblas-base  0.1alpha2.2-3          Optimized BLAS (linear algebra) library based on GotoBLAS2 
!ii  libopenblas-dev   0.1alpha2.2-3          Optimized BLAS (linear algebra) library based on GotoBLAS2 
!>
 

Ada

-- Ada implementation 
with Ada.Text_Io; use Ada.Text_Io; 
with Ada.Float_Text_Io; use Ada.Float_Text_Io; 
with Ada.integer_Text_Io; use Ada.integer_Text_Io; 
with Ada.Numerics.Real_Arrays;  use Ada.Numerics.Real_Arrays; 
procedure t44 is 
    A : constant real_matrix := 
             (( 1.0,  2.0,  3.0), 
              ( 4.0,  5.0,  6.0), 
              ( 7.0,  8.0,  9.0)); 
           --(( -3.0,  2.0,  -5.0), 
           --   ( -1.0,  0.0,  -2.0), 
           --   ( 3.0,  -4.0,  1.0)); 
    C : real_matrix(A'range(1), A'range(2)); 
    B : real_matrix(1 .. 2, 1 .. 2); 
    -- Thanks goes to Dmitry A. Kazakov this function 
    function Exclude (A : Real_Matrix; I, J : Integer) 
       return Real_Matrix is 
       AI, AJ : Integer := A'First (1); 
    begin 
       return B : Real_Matrix (1..A'Length (1)-1, 1..A'Length (2)-1) do 
          AI := A'First (1); 
          for BI in B'Range (1) loop 
             if AI = I then 
                AI := AI + 1; 
             end if; 
             AJ := A'First (2); 
             for BJ in B'Range (2) loop 
                if AJ = J then 
                   AJ := AJ + 1; 
                end if; 
                B (BI, BJ) := A (AI, AJ); 
                AJ := AJ + 1; 
             end loop; 
             AI := AI + 1; 
          end loop; 
       end return; 
    end Exclude; 
    -- Function to print a 2D matrix 
    procedure Put (X : Real_Matrix) is 
       begin 
          for I in X'Range (1) loop 
             for J in X'Range (2) loop 
                Put (X (I, J) ); put("    "); 
             end loop; 
             New_Line; 
          end loop; 
   end Put; 
begin 
    FOR I in A'range(1) LOOP 
      FOR J in A'range(2) LOOP 
          B := Exclude (A, I, J); 
          C(I,J) := float((-1)**( ((I+J) rem 2) ))*determinant(B); 
      END LOOP; 
    END LOOP; 
    c := transpose(c); 
    put(c); 
end t44;

compile and run

-- export LD_LIBRARY_PATH=/usr/lib/atlas-base/:/usr/lib/atlas-base/atlas/ 
>gnatmake t44.adb 
gcc-4.6 -c t44.adb 
gnatbind -x t44.ali 
gnatlink t44.ali 
 
>./t44 
-3.00000E+00     6.00000E+00    -3.00000E+00 
 6.00000E+00    -1.20000E+01     6.00000E+00 
-3.00000E+00     6.00000E+00    -3.00000E+00