Mercurial > hg > octave-lyh
annotate src/DLD-FUNCTIONS/kron.cc @ 10527:b4d2080b6df7
Replace nzmax by nnz as needed
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Fri, 16 Apr 2010 16:14:45 +0200 |
| parents | 81067c72361f |
| children | 89f4d7e294cc |
| rev | line source |
|---|---|
| 3910 | 1 /* |
| 2 | |
| 8920 | 3 Copyright (C) 2002, 2005, 2006, 2007, 2008 John W. Eaton |
| 3910 | 4 |
| 5 This file is part of Octave. | |
| 6 | |
| 7 Octave is free software; you can redistribute it and/or modify it | |
| 8 under the terms of the GNU General Public License as published by the | |
| 7016 | 9 Free Software Foundation; either version 3 of the License, or (at your |
| 10 option) any later version. | |
| 3910 | 11 |
| 12 Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 15 for more details. | |
| 16 | |
| 17 You should have received a copy of the GNU General Public License | |
| 7016 | 18 along with Octave; see the file COPYING. If not, see |
| 19 <http://www.gnu.org/licenses/>. | |
| 3910 | 20 |
| 21 */ | |
| 22 | |
| 3911 | 23 // Author: Paul Kienzle <pkienzle@users.sf.net> |
| 24 | |
| 3910 | 25 #ifdef HAVE_CONFIG_H |
| 26 #include <config.h> | |
| 27 #endif | |
| 28 | |
| 29 #include "dMatrix.h" | |
| 10461 | 30 #include "fMatrix.h" |
| 3910 | 31 #include "CMatrix.h" |
| 10461 | 32 #include "fCMatrix.h" |
| 33 | |
| 34 #include "dSparse.h" | |
| 35 #include "CSparse.h" | |
| 36 | |
| 37 #include "dDiagMatrix.h" | |
| 38 #include "fDiagMatrix.h" | |
| 39 #include "CDiagMatrix.h" | |
| 40 #include "fCDiagMatrix.h" | |
| 41 | |
| 42 #include "PermMatrix.h" | |
| 43 | |
| 44 #include "mx-inlines.cc" | |
| 4153 | 45 #include "quit.h" |
| 3910 | 46 |
| 47 #include "defun-dld.h" | |
| 48 #include "error.h" | |
| 49 #include "oct-obj.h" | |
| 50 | |
| 10461 | 51 template <class R, class T> |
| 52 static MArray<T> | |
| 53 kron (const MArray<R>& a, const MArray<T>& b) | |
| 54 { | |
| 55 assert (a.ndims () == 2); | |
| 56 assert (b.ndims () == 2); | |
| 57 | |
| 58 octave_idx_type nra = a.rows (), nrb = b.rows (); | |
| 59 octave_idx_type nca = a.cols (), ncb = b.cols (); | |
| 3910 | 60 |
| 10461 | 61 MArray<T> c (nra*nrb, nca*ncb); |
| 62 T *cv = c.fortran_vec (); | |
| 63 | |
| 64 for (octave_idx_type ja = 0; ja < nca; ja++) | |
| 65 for (octave_idx_type jb = 0; jb < ncb; jb++) | |
| 66 for (octave_idx_type ia = 0; ia < nra; ia++) | |
| 67 { | |
| 68 octave_quit (); | |
| 69 mx_inline_mul (nrb, cv, a(ia, ja), b.data () + nrb*jb); | |
| 70 cv += nrb; | |
| 71 } | |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7455
diff
changeset
|
72 |
| 10461 | 73 return c; |
| 74 } | |
| 75 | |
| 76 template <class R, class T> | |
| 77 static MArray<T> | |
| 78 kron (const MDiagArray2<R>& a, const MArray<T>& b) | |
| 79 { | |
| 80 assert (b.ndims () == 2); | |
| 81 | |
| 82 octave_idx_type nra = a.rows (), nrb = b.rows (), dla = a.diag_length (); | |
| 83 octave_idx_type nca = a.cols (), ncb = b.cols (); | |
|
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7455
diff
changeset
|
84 |
| 10461 | 85 MArray<T> c (nra*nrb, nca*ncb, T()); |
| 86 | |
| 87 for (octave_idx_type ja = 0; ja < dla; ja++) | |
| 88 for (octave_idx_type jb = 0; jb < ncb; jb++) | |
| 89 { | |
| 90 octave_quit (); | |
| 91 mx_inline_mul (nrb, &c.xelem(ja*nrb, ja*ncb + jb), a.dgelem (ja), b.data () + nrb*jb); | |
| 92 } | |
| 93 | |
| 94 return c; | |
| 95 } | |
| 3910 | 96 |
| 97 template <class T> | |
| 10461 | 98 static MSparse<T> |
| 99 kron (const MSparse<T>& A, const MSparse<T>& B) | |
|
7455
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
100 { |
|
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
101 octave_idx_type idx = 0; |
| 10461 | 102 MSparse<T> C (A.rows () * B.rows (), A.columns () * B.columns (), |
|
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10461
diff
changeset
|
103 A.nnz () * B.nnz ()); |
|
7455
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
104 |
|
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
105 C.cidx (0) = 0; |
|
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
106 |
|
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
107 for (octave_idx_type Aj = 0; Aj < A.columns (); Aj++) |
|
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
108 for (octave_idx_type Bj = 0; Bj < B.columns (); Bj++) |
|
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
109 { |
| 10461 | 110 octave_quit (); |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
111 for (octave_idx_type Ai = A.cidx (Aj); Ai < A.cidx (Aj+1); Ai++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
112 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
113 octave_idx_type Ci = A.ridx(Ai) * B.rows (); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
114 const T v = A.data (Ai); |
|
7455
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
115 |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
116 for (octave_idx_type Bi = B.cidx (Bj); Bi < B.cidx (Bj+1); Bi++) |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
117 { |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
118 C.data (idx) = v * B.data (Bi); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
119 C.ridx (idx++) = Ci + B.ridx (Bi); |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
120 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
121 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
122 C.cidx (Aj * B.columns () + Bj + 1) = idx; |
|
7455
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
123 } |
| 10461 | 124 |
| 125 return C; | |
|
7455
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
126 } |
|
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
127 |
| 10461 | 128 static PermMatrix |
| 129 kron (const PermMatrix& a, const PermMatrix& b) | |
| 130 { | |
| 131 octave_idx_type na = a.rows (), nb = b.rows (); | |
| 132 const octave_idx_type *pa = a.data (), *pb = b.data (); | |
| 133 PermMatrix c(na*nb); // Row permutation. | |
| 134 octave_idx_type *pc = c.fortran_vec (); | |
|
7455
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
135 |
| 10461 | 136 bool cola = a.is_col_perm (), colb = b.is_col_perm (); |
| 137 if (cola && colb) | |
| 138 { | |
| 139 for (octave_idx_type i = 0; i < na; i++) | |
| 140 for (octave_idx_type j = 0; j < nb; j++) | |
| 141 pc[pa[i]*nb+pb[j]] = i*nb+j; | |
| 142 } | |
| 143 else if (cola) | |
| 144 { | |
| 145 for (octave_idx_type i = 0; i < na; i++) | |
| 146 for (octave_idx_type j = 0; j < nb; j++) | |
| 147 pc[pa[i]*nb+j] = i*nb+pb[j]; | |
| 148 } | |
| 149 else if (colb) | |
| 150 { | |
| 151 for (octave_idx_type i = 0; i < na; i++) | |
| 152 for (octave_idx_type j = 0; j < nb; j++) | |
| 153 pc[i*nb+pb[j]] = pa[i]*nb+j; | |
| 154 } | |
| 155 else | |
| 156 { | |
| 157 for (octave_idx_type i = 0; i < na; i++) | |
| 158 for (octave_idx_type j = 0; j < nb; j++) | |
| 159 pc[i*nb+j] = pa[i]*nb+pb[j]; | |
| 160 } | |
| 161 | |
| 162 return c; | |
| 163 } | |
|
7455
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
164 |
|
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
165 |
| 10461 | 166 template <class MTA, class MTB> |
| 167 octave_value | |
| 168 do_kron (const octave_value& a, const octave_value& b) | |
| 169 { | |
| 170 MTA am = octave_value_extract<MTA> (a); | |
| 171 MTB bm = octave_value_extract<MTB> (b); | |
| 172 return octave_value (kron (am, bm)); | |
| 173 } | |
| 174 | |
| 175 #define ALL_TYPES(AMT, BMT) \ | |
| 176 } while (0) \ | |
| 177 | |
| 178 DEFUN_DLD (kron, args, , "-*- texinfo -*-\n\ | |
| 6678 | 179 @deftypefn {Loadable Function} {} kron (@var{a}, @var{b})\n\ |
| 3910 | 180 Form the kronecker product of two matrices, defined block by block as\n\ |
| 181 \n\ | |
| 182 @example\n\ | |
| 183 x = [a(i, j) b]\n\ | |
| 184 @end example\n\ | |
| 185 \n\ | |
| 186 For example,\n\ | |
| 187 \n\ | |
| 188 @example\n\ | |
| 189 @group\n\ | |
| 190 kron (1:4, ones (3, 1))\n\ | |
| 191 @result{} 1 2 3 4\n\ | |
| 192 1 2 3 4\n\ | |
| 193 1 2 3 4\n\ | |
| 194 @end group\n\ | |
| 195 @end example\n\ | |
| 196 @end deftypefn") | |
| 197 { | |
| 10461 | 198 octave_value retval; |
| 3910 | 199 |
| 200 int nargin = args.length (); | |
| 201 | |
| 10461 | 202 if (nargin == 2) |
| 3910 | 203 { |
| 10461 | 204 octave_value a = args(0), b = args(1); |
| 205 if (a.is_perm_matrix () && b.is_perm_matrix ()) | |
| 206 retval = do_kron<PermMatrix, PermMatrix> (a, b); | |
| 207 else if (a.is_diag_matrix ()) | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
208 { |
| 10461 | 209 if (b.is_diag_matrix () && a.rows () == a.columns () |
| 210 && b.rows () == b.columns ()) | |
| 211 { | |
| 212 octave_value_list tmp; | |
| 213 tmp(0) = a.diag (); | |
| 214 tmp(1) = b.diag (); | |
| 215 tmp = Fkron (tmp, 1); | |
| 216 if (tmp.length () == 1) | |
| 217 retval = tmp(0).diag (); | |
| 218 } | |
| 219 else if (a.is_single_type () || b.is_single_type ()) | |
| 220 { | |
| 221 if (a.is_complex_type ()) | |
| 222 return do_kron<FloatComplexDiagMatrix, FloatComplexMatrix> (a, b); | |
| 223 else if (b.is_complex_type ()) | |
| 224 return do_kron<FloatDiagMatrix, FloatComplexMatrix> (a, b); | |
| 225 else | |
| 226 return do_kron<FloatDiagMatrix, FloatMatrix> (a, b); | |
| 227 } | |
| 228 else | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
229 { |
| 10461 | 230 if (a.is_complex_type ()) |
| 231 return do_kron<ComplexDiagMatrix, ComplexMatrix> (a, b); | |
| 232 else if (b.is_complex_type ()) | |
| 233 return do_kron<DiagMatrix, ComplexMatrix> (a, b); | |
| 234 else | |
| 235 return do_kron<DiagMatrix, Matrix> (a, b); | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
236 } |
|
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
237 } |
| 10461 | 238 else if (a.is_sparse_type () || b.is_sparse_type ()) |
| 239 { | |
| 240 if (args(0).is_complex_type () || args(1).is_complex_type ()) | |
| 241 return do_kron<SparseComplexMatrix, SparseComplexMatrix> (a, b); | |
| 242 else | |
| 243 return do_kron<SparseMatrix, SparseMatrix> (a, b); | |
| 244 } | |
| 245 else if (a.is_single_type () || b.is_single_type ()) | |
| 246 { | |
| 247 if (a.is_complex_type ()) | |
| 248 return do_kron<FloatComplexMatrix, FloatComplexMatrix> (a, b); | |
| 249 else if (b.is_complex_type ()) | |
| 250 return do_kron<FloatMatrix, FloatComplexMatrix> (a, b); | |
| 251 else | |
| 252 return do_kron<FloatMatrix, FloatMatrix> (a, b); | |
| 253 } | |
|
7455
fe332ce262b5
eliminate spkron.cc; dispatch in kron
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
254 else |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
255 { |
| 10461 | 256 if (a.is_complex_type ()) |
| 257 return do_kron<ComplexMatrix, ComplexMatrix> (a, b); | |
| 258 else if (b.is_complex_type ()) | |
| 259 return do_kron<Matrix, ComplexMatrix> (a, b); | |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
260 else |
| 10461 | 261 return do_kron<Matrix, Matrix> (a, b); |
|
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9209
diff
changeset
|
262 } |
| 3910 | 263 } |
| 264 | |
| 265 return retval; | |
| 266 } |