### 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}$

ﬁnd 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 ﬁnd function in Matlab to do this. But for non-singular matrices only the direct method of ﬁnding 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 implementation
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/
`