octave-nkf: liboctave/oct-convn.cc annotate

author	Jaroslav Hajek <highegg@gmail.com>
date	Tue, 02 Mar 2010 10:59:05 +0100
parents
children	5af0b4bb384d

rev	line source
10385 56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	1 /*
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	2
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	3 Copyright (C) 2010 VZLU Prague
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	4
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	5 This file is part of Octave.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	6
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	7 Octave is free software; you can redistribute it and/or modify it
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	8 under the terms of the GNU General Public License as published by the
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	9 Free Software Foundation; either version 3 of the License, or (at your
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	10 option) any later version.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	11
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	12 Octave is distributed in the hope that it will be useful, but WITHOUT
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	15 for more details.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	16
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	17 You should have received a copy of the GNU General Public License
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	18 along with Octave; see the file COPYING. If not, see
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	19 <http://www.gnu.org/licenses/>.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	20
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	21 */
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	22
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	23 #ifdef HAVE_CONFIG_H
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	24 #include <config.h>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	25 #endif
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	26
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	27 #include <iostream>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	28 #include <algorithm>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	29 #include "oct-convn.h"
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	30 #include "oct-locbuf.h"
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	31
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	32
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	33 // FIXME: Only the axpy-form accumulation is used. This is natural for outer
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	34 // convolution as it requires no boundary treatment.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	35 // This typically requires one more memory store per operation, but as the
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	36 // memory access pattern is equivalent (switching ro and rw parts), I wouldn't
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	37 // expect a significant difference. cf. for instance sum(), where row-wise sum
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	38 // (axpy-based) is no slower than column-wise (dot-based).
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	39 // It would be nice, however, if the inner convolution was computed directly by
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	40 // dot-based accumulation.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	41
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	42 // FIXME: Specifying the kernel as outer product should probably get special
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	43 // treatment.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	44
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	45 // All kernels smaller than 7x5 get specialized code.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	46 #define MAX_KERNEL_SIZE_M 7
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	47 #define MAX_KERNEL_SIZE_N 5
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	48
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	49 template <class T, class R, int mb, int nb>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	50 static void
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	51 convolve_2d_kernel_axpy (const T *a, octave_idx_type ma, octave_idx_type na,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	52 const R b, T c, octave_idx_type ldc)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	53 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	54 for (octave_idx_type ja = 0; ja < na; ja++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	55 for (octave_idx_type ia = 0; ia < ma; ia++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	56 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	57 T aij = a[ma*ja + ia];
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	58 // The following double loop is a candidate for complete unrolling.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	59 for (int jb = 0; jb < nb; jb++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	60 for (int ib = 0; ib < mb; ib++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	61 c[(ja+jb)ldc + (ia+ib)] += aij b[mb*jb+ib];
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	62 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	63 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	64
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	65 // Kernel dispatcher.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	66 template <class T, class R>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	67 static void
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	68 convolve_2d_axpy (const T *a, octave_idx_type ma, octave_idx_type na,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	69 const R *b, octave_idx_type mb, octave_idx_type nb,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	70 T *c, octave_idx_type ldc)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	71 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	72 // Table of kernels.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	73 static void (*table[MAX_KERNEL_SIZE_M][MAX_KERNEL_SIZE_N])
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	74 (const T , octave_idx_type, octave_idx_type, const R , T *, octave_idx_type)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	75 = {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	76 // This must be repeated MAX_KERNEL_SIZE-times. I do not see a way to
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	77 // automate this.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	78 #define STORE_KERNEL_POINTER(M,N) \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	79 convolve_2d_kernel_axpy<T,R,M,N>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	80 #define STORE_KERNEL_ROW(M) \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	81 { \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	82 STORE_KERNEL_POINTER(M,1), \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	83 STORE_KERNEL_POINTER(M,2), \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	84 STORE_KERNEL_POINTER(M,3), \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	85 STORE_KERNEL_POINTER(M,4), \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	86 STORE_KERNEL_POINTER(M,5), \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	87 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	88
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	89 STORE_KERNEL_ROW(1),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	90 STORE_KERNEL_ROW(2),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	91 STORE_KERNEL_ROW(3),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	92 STORE_KERNEL_ROW(4),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	93 STORE_KERNEL_ROW(5),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	94 STORE_KERNEL_ROW(6),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	95 STORE_KERNEL_ROW(7),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	96 };
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	97
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	98 if (mb != 0 && nb != 0)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	99 (*table[mb-1][nb-1]) (a, ma, na, b, c, ldc);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	100 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	101
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	102 // 2d convolution with a matrix kernel.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	103 template <class T, class R>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	104 static void
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	105 convolve_2d (const T *a, octave_idx_type ma, octave_idx_type na,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	106 const R *b, octave_idx_type mb, octave_idx_type nb,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	107 T *c)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	108 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	109 octave_idx_type ldc = ma + mb - 1;
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	110 if (mb <= MAX_KERNEL_SIZE_M && nb <= MAX_KERNEL_SIZE_N)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	111 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	112 // Call kernel directly on b.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	113 convolve_2d_axpy (a, ma, na, b, mb, nb, c, ldc);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	114 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	115 else
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	116 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	117 // Split b to blocks.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	118 OCTAVE_LOCAL_BUFFER (R, b1, MAX_KERNEL_SIZE_M*MAX_KERNEL_SIZE_N);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	119 for (octave_idx_type jb = 0; jb < nb; jb += MAX_KERNEL_SIZE_N)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	120 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	121 octave_idx_type nb1 = std::min (nb - jb, MAX_KERNEL_SIZE_N);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	122 for (octave_idx_type ib = 0; ib < mb; ib += MAX_KERNEL_SIZE_M)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	123 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	124 octave_idx_type mb1 = std::min (mb - ib, MAX_KERNEL_SIZE_M);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	125
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	126 // Copy block to buffer.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	127 R *bf = b1;
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	128 for (octave_idx_type j = jb; j < jb+nb1; j++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	129 for (octave_idx_type i = ib; i < ib+mb1; i++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	130 bf++ = b[jmb + i];
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	131
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	132 // Call kernel.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	133 convolve_2d_axpy (a, ma, na, b1, mb1, nb1, c + ldc*jb + ib, ldc);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	134 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	135 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	136 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	137
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	138 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	139
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	140 template <class T, class R>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	141 void convolve_nd (const T *a, const dim_vector& ad, const dim_vector& acd,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	142 const R *b, const dim_vector& bd, const dim_vector& bcd,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	143 T *c, const dim_vector& ccd, int nd)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	144 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	145 if (nd == 2)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	146 convolve_2d<T, R> (a, ad(0), ad(1), b, bd(0), bd(1), c);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	147 else
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	148 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	149 octave_idx_type ma = acd(nd-2), na = ad(nd-1), mb = bcd(nd-2), nb = bd(nd-1);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	150 octave_idx_type ldc = ccd(nd-2);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	151 for (octave_idx_type jb = 0; jb < nb; jb++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	152 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	153 for (octave_idx_type ja = 0; ja < na; ja++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	154 convolve_nd<T, R> (a + maja, ad, acd, b + mbjb, bd, bcd,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	155 c + ldc*(ja+jb), ccd, nd-1);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	156 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	157 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	158 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	159
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	160 // Arbitrary convolutor.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	161 // The 2nd array is assumed to be the smaller one.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	162 template <class T, class R>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	163 static MArray<T>
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	164 convolve (const MArray<T>& a, const MArray<R>& b,
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	165 convn_type ct)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	166 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	167 if (a.is_empty () \|\| b.is_empty ())
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	168 return MArray<T> ();
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	169
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	170 int nd = std::max (a.ndims (), b.ndims ());
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	171 const dim_vector adims = a.dims ().redim (nd), bdims = b.dims ().redim (nd);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	172 dim_vector cdims = dim_vector::alloc (nd);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	173
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	174 for (int i = 0; i < nd; i++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	175 cdims(i) = std::max (adims(i) + bdims(i) - 1, 0);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	176
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	177 MArray<T> c (cdims, T());
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	178
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	179 convolve_nd<T, R> (a.fortran_vec (), adims, adims.cumulative (),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	180 b.fortran_vec (), bdims, bdims.cumulative (),
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	181 c.fortran_vec (), cdims.cumulative (), nd);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	182
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	183 // Pick the relevant part.
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	184 Array<idx_vector> sidx (nd, 1);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	185
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	186 switch (ct)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	187 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	188 case convn_valid:
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	189 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	190 for (int i = 0; i < nd; i++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	191 sidx(i) = idx_vector (bdims(i)-1, adims(i));
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	192 c = c.index (sidx);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	193 break;
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	194 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	195 case convn_same:
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	196 {
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	197 for (int i = 0; i < nd; i++)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	198 sidx(i) = idx_vector::make_range ((bdims(i)-1)/2, 1, adims(i));
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	199 c = c.index (sidx);
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	200 break;
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	201 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	202 default:
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	203 break;
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	204 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	205
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	206 return c;
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	207 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	208
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	209 #define CONV_DEFS(TPREF, RPREF) \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	210 TPREF ## NDArray \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	211 convn (const TPREF ## NDArray& a, const RPREF ## NDArray& b, convn_type ct) \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	212 { \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	213 return convolve (a, b, ct); \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	214 } \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	215 TPREF ## Matrix \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	216 convn (const TPREF ## Matrix& a, const RPREF ## Matrix& b, convn_type ct) \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	217 { \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	218 return convolve (a, b, ct); \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	219 } \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	220 TPREF ## Matrix \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	221 convn (const TPREF ## Matrix& a, const RPREF ## ColumnVector& c, \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	222 const RPREF ## RowVector& r, convn_type ct) \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	223 { \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	224 return convolve (a, c * r, ct); \
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	225 }
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	226
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	227 CONV_DEFS ( , )
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	228 CONV_DEFS (Complex, )
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	229 CONV_DEFS (Complex, Complex)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	230 CONV_DEFS (Float, Float)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	231 CONV_DEFS (FloatComplex, Float)
56116dceb1e0 add omitted source from the last change Jaroslav Hajek <highegg@gmail.com> parents: diff changeset	232 CONV_DEFS (FloatComplex, FloatComplex)

10385

56116dceb1e0 add omitted source from the last change