2540
|
1 ## Copyright (C) 1995, 1996 Kurt Hornik |
3426
|
2 ## |
3922
|
3 ## This file is part of Octave. |
|
4 ## |
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
6 ## under the terms of the GNU General Public License as published by |
7016
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
8 ## your option) any later version. |
3426
|
9 ## |
3922
|
10 ## Octave is distributed in the hope that it will be useful, but |
2540
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
3426
|
13 ## General Public License for more details. |
|
14 ## |
2540
|
15 ## You should have received a copy of the GNU General Public License |
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
17 ## <http://www.gnu.org/licenses/>. |
2540
|
18 |
3321
|
19 ## -*- texinfo -*- |
|
20 ## @deftypefn {Function File} {} commutation_matrix (@var{m}, @var{n}) |
|
21 ## Return the commutation matrix |
|
22 ## @iftex |
|
23 ## @tex |
|
24 ## $K_{m,n}$ |
|
25 ## @end tex |
|
26 ## @end iftex |
|
27 ## @ifinfo |
|
28 ## K(m,n) |
|
29 ## @end ifinfo |
|
30 ## which is the unique |
|
31 ## @iftex |
|
32 ## @tex |
|
33 ## $m n \times m n$ |
|
34 ## @end tex |
|
35 ## @end iftex |
|
36 ## @ifinfo |
|
37 ## @var{m}*@var{n} by @var{m}*@var{n} |
|
38 ## @end ifinfo |
|
39 ## matrix such that |
|
40 ## @iftex |
|
41 ## @tex |
|
42 ## $K_{m,n} \cdot {\rm vec} (A) = {\rm vec} (A^T)$ |
|
43 ## @end tex |
|
44 ## @end iftex |
|
45 ## @ifinfo |
3499
|
46 ## @math{K(m,n) * vec(A) = vec(A')} |
3321
|
47 ## @end ifinfo |
|
48 ## for all |
|
49 ## @iftex |
|
50 ## @tex |
|
51 ## $m\times n$ |
|
52 ## @end tex |
|
53 ## @end iftex |
|
54 ## @ifinfo |
3499
|
55 ## @math{m} by @math{n} |
3321
|
56 ## @end ifinfo |
|
57 ## matrices |
|
58 ## @iftex |
|
59 ## @tex |
|
60 ## $A$. |
|
61 ## @end tex |
|
62 ## @end iftex |
|
63 ## @ifinfo |
3499
|
64 ## @math{A}. |
3321
|
65 ## @end ifinfo |
3426
|
66 ## |
3321
|
67 ## If only one argument @var{m} is given, |
|
68 ## @iftex |
|
69 ## @tex |
|
70 ## $K_{m,m}$ |
|
71 ## @end tex |
|
72 ## @end iftex |
|
73 ## @ifinfo |
3499
|
74 ## @math{K(m,m)} |
3321
|
75 ## @end ifinfo |
|
76 ## is returned. |
3426
|
77 ## |
2540
|
78 ## See Magnus and Neudecker (1988), Matrix differential calculus with |
|
79 ## applications in statistics and econometrics. |
3321
|
80 ## @end deftypefn |
2540
|
81 |
5428
|
82 ## Author: KH <Kurt.Hornik@wu-wien.ac.at> |
2540
|
83 ## Created: 8 May 1995 |
|
84 ## Adapted-By: jwe |
|
85 |
|
86 function k = commutation_matrix (m, n) |
3426
|
87 |
2540
|
88 if (nargin < 1 || nargin > 2) |
6046
|
89 print_usage (); |
2540
|
90 else |
4030
|
91 if (! (isscalar (m) && m == round (m) && m > 0)) |
2540
|
92 error ("commutation_matrix: m must be a positive integer"); |
|
93 endif |
|
94 if (nargin == 1) |
|
95 n = m; |
4030
|
96 elseif (! (isscalar (n) && n == round (n) && n > 0)) |
2540
|
97 error ("commutation_matrix: n must be a positive integer"); |
|
98 endif |
|
99 endif |
3426
|
100 |
2540
|
101 ## It is clearly possible to make this a LOT faster! |
|
102 k = zeros (m * n, m * n); |
|
103 for i = 1 : m |
|
104 for j = 1 : n |
|
105 k ((i - 1) * n + j, (j - 1) * m + i) = 1; |
|
106 endfor |
|
107 endfor |
|
108 |
|
109 endfunction |