diff --git a/py/fastspecfit/continuum.py b/py/fastspecfit/continuum.py
index 6f5ac188..9c24e714 100644
--- a/py/fastspecfit/continuum.py
+++ b/py/fastspecfit/continuum.py
@@ -539,7 +539,7 @@ def build_stellar_continuum(self, templateflux, templatecoeff,
 
         Parameters
         ----------
-        templateflux : :class:`numpy.ndarray` [npix, ntemplates]
+        templateflux : :class:`numpy.ndarray` [ntemplates, npix]
             Rest-frame, native-resolution template spectra corresponding to
             `templatewave`.
         templatecoeff : :class:`numpy.ndarray` [ntemplates]
@@ -566,7 +566,7 @@ def build_stellar_continuum(self, templateflux, templatecoeff,
         """
         if conv_pre is None or vdisp > Templates.MAX_PRE_VDISP:
             # Compute the weighted sum of the templates.
-            contmodel = templateflux.dot(templatecoeff)
+            contmodel = templatecoeff.dot(templateflux)
 
             # Optionally convolve to the desired velocity dispersion.
             if vdisp is not None:
@@ -578,13 +578,13 @@ def build_stellar_continuum(self, templateflux, templatecoeff,
             flux_lohi, ft_flux_mid, fft_len = conv_pre
 
             # Compute the weighted sum of the templates.
-            cont_lohi   = flux_lohi.dot(templatecoeff)
-            ft_cont_mid = ft_flux_mid.dot(templatecoeff)
+            cont_lohi   = templatecoeff.dot(flux_lohi)
+            ft_cont_mid = templatecoeff.dot(ft_flux_mid)
 
             # Convolve to the desired velocity dispersion. Use the vdisp
             # convolution that takes precomputed FT of flux for convolved
             # region.
-            flux_len = templateflux.shape[0]
+            flux_len = templateflux.shape[1]
             contmodel = self.templates.convolve_vdisp_from_pre(
                 cont_lohi, ft_cont_mid, flux_len, fft_len, vdisp)
 
@@ -759,7 +759,7 @@ def fit_stellar_continuum(self, templateflux, fit_vdisp=False, conv_pre=None,
 
         Parameters
         ----------
-        templateflux : :class:`numpy.ndarray` [npix, ntemplate]
+        templateflux : :class:`numpy.ndarray` [ntemplate, npix]
             Grid of input (model) spectra.
         fit_vdisp : :class:`bool`
             If `True`, solve for the velocity dispersion;
@@ -830,7 +830,7 @@ def fit_stellar_continuum(self, templateflux, fit_vdisp=False, conv_pre=None,
         if vdisp_bounds is None:
             vdisp_bounds = self.vdisp_bounds
 
-        ntemplates = templateflux.shape[1]
+        ntemplates = templateflux.shape[0]
 
         # Unpack the input data to infer the fitting "mode" and the objective
         # function.
@@ -995,7 +995,7 @@ def continuum_fastphot(redshift, objflam, objflamivar, CTools,
 
         t0 = time.time()
         tauv, _, coeff, resid = CTools.fit_stellar_continuum(
-            templates.flux_nomvdisp[:, agekeep], fit_vdisp=False,
+            templates.flux_nomvdisp[agekeep, :], fit_vdisp=False,
             objflam=objflam, objflamistd=objflamistd,
             synthphot=True, synthspec=False)
         dt = time.time()-t0
@@ -1017,7 +1017,7 @@ def continuum_fastphot(redshift, objflam, objflamivar, CTools,
             # Get the best-fitting model with and without line-emission.
             sedmodel = CTools.optimizer_saved_contmodel
             sedmodel_nolines = CTools.build_stellar_continuum(
-                templates.flux_nolines_nomvdisp[:, agekeep], coeff,
+                templates.flux_nolines_nomvdisp[agekeep, :], coeff,
                 tauv=tauv, vdisp=None, dust_emission=False)
 
             # Measure Dn(4000) from the line-free model.
@@ -1042,7 +1042,7 @@ def continuum_fastphot(redshift, objflam, objflamivar, CTools,
 
             for imonte in range(nmonte):
                 tauv1, _, coeff1, _ = CTools.fit_stellar_continuum(
-                    templates.flux_nomvdisp[:, agekeep], fit_vdisp=False,
+                    templates.flux_nomvdisp[agekeep, :], fit_vdisp=False,
                     tauv_guess=tauv_guess[imonte], objflam=objflam_monte[:, imonte],
                     objflamistd=objflamistd, synthphot=True, synthspec=False)
 
@@ -1051,7 +1051,7 @@ def continuum_fastphot(redshift, objflam, objflamivar, CTools,
 
                 sedmodel_monte[:, imonte] = CTools.optimizer_saved_contmodel
                 sedmodel_nolines_monte[:, imonte] = CTools.build_stellar_continuum(
-                    templates.flux_nolines_nomvdisp[:, agekeep], coeff1,
+                    templates.flux_nolines_nomvdisp[agekeep, :], coeff1,
                     tauv=tauv1, vdisp=None, dust_emission=False)
 
                 dn4000_model1, _ = Photometry.get_dn4000(
@@ -1158,7 +1158,7 @@ def continuum_fastspec(redshift, objflam, objflamivar, CTools,
     # the IR fluxes put an additional constraint on the dust content.
     init_tauv_bounds = (0., 1.)
     tauv_nomvdisp, _, coeff_nomvdisp, resid_nomvdisp = CTools.fit_stellar_continuum(
-        templates.flux_nolines_nomvdisp[:, agekeep], fit_vdisp=False, conv_pre=None,
+        templates.flux_nolines_nomvdisp[agekeep, :], fit_vdisp=False, conv_pre=None,
         tauv_bounds=init_tauv_bounds, specflux=specflux, specistd=specistd,
         dust_emission=False, synthspec=True)
 
@@ -1178,7 +1178,7 @@ def continuum_fastspec(redshift, objflam, objflamivar, CTools,
         input_conv_pre_nolines = templates.conv_pre_select(templates.conv_pre_nolines, agekeep)
 
         tauv_withvdisp, vdisp, coeff_withvdisp, resid_withvdisp = CTools.fit_stellar_continuum(
-            templates.flux_nolines[:, agekeep], fit_vdisp=True,
+            templates.flux_nolines[agekeep, :], fit_vdisp=True,
             conv_pre=input_conv_pre_nolines,
             vdisp_guess=templates.vdisp_nominal,
             tauv_bounds=init_tauv_bounds, specflux=specflux,
@@ -1198,8 +1198,8 @@ def continuum_fastspec(redshift, objflam, objflamivar, CTools,
             contmodel = contmodel_withvdisp
 
             # convolve the templates at the derived vdisp
-            input_templateflux = templates.convolve_vdisp(templates.flux[:, agekeep], vdisp)
-            input_templateflux_nolines = templates.convolve_vdisp(templates.flux_nolines[:, agekeep], vdisp)
+            input_templateflux = templates.convolve_vdisp(templates.flux[agekeep, :], vdisp)
+            input_templateflux_nolines = templates.convolve_vdisp(templates.flux_nolines[agekeep, :], vdisp)
         else:
             tauv = tauv_nomvdisp
             coeff = coeff_nomvdisp
@@ -1219,7 +1219,7 @@ def continuum_fastspec(redshift, objflam, objflamivar, CTools,
 
             for imonte in range(nmonte):
                 tauv1, vdisp1, coeff1, _ = CTools.fit_stellar_continuum(
-                    templates.flux_nolines[:, agekeep], fit_vdisp=True,
+                    templates.flux_nolines[agekeep, :], fit_vdisp=True,
                     conv_pre=input_conv_pre_nolines,
                     vdisp_guess=vdisp_guess[imonte],
                     tauv_guess=tauv_guess[imonte],
@@ -1302,8 +1302,8 @@ def continuum_fastspec(redshift, objflam, objflamivar, CTools,
         tauv = tauv_nomvdisp
         coeff = coeff_nomvdisp
         vdisp = templates.vdisp_nominal
-        input_templateflux = templates.flux_nomvdisp[:, agekeep]
-        input_templateflux_nolines = templates.flux_nolines_nomvdisp[:, agekeep]
+        input_templateflux = templates.flux_nomvdisp[agekeep, :]
+        input_templateflux_nolines = templates.flux_nolines_nomvdisp[agekeep, :]
         contmodel = CTools.optimizer_saved_contmodel.copy()
 
     # Next, estimate the aperture correction.
diff --git a/py/fastspecfit/templates.py b/py/fastspecfit/templates.py
index 35fa78e7..cb95c3fe 100644
--- a/py/fastspecfit/templates.py
+++ b/py/fastspecfit/templates.py
@@ -76,6 +76,9 @@ def __init__(self, template_file=None, template_version=None, imf=None,
         templateinfo     = T['METADATA'].read()
         templatehdr      = T['METADATA'].read_header()
 
+        templateflux     = np.transpose(templateflux).copy()
+        templatelineflux = np.transpose(templatelineflux).copy()
+
         self.version = T[0].read_header()['VERSION']
 
         self.imf = templatehdr['IMF']
@@ -87,8 +90,8 @@ def __init__(self, template_file=None, template_version=None, imf=None,
         keeplo = np.searchsorted(templatewave, mintemplatewave, 'left')
         keephi = np.searchsorted(templatewave, maxtemplatewave, 'right')
         self.wave = templatewave[keeplo:keephi]
-        self.flux = templateflux[keeplo:keephi, :]
-        self.flux_nolines = self.flux - templatelineflux[keeplo:keephi, :]
+        self.flux = templateflux[:, keeplo:keephi]
+        self.flux_nolines = self.flux - templatelineflux[:, keeplo:keephi]
         self.npix = len(self.wave)
 
         # dust attenuation curve
@@ -196,7 +199,7 @@ def convolve_vdisp_pre(self, templateflux):
         Parameters
         ----------
         templateflux: :class:`np.ndarray`
-            [nwavelengths x ntemplates] array of templates
+            [ntemplates x nwavelengths] array of templates
 
         Returns
         -------
@@ -215,43 +218,43 @@ def convolve_vdisp_pre(self, templateflux):
         lo, hi = self.pixkms_bounds
 
         # extract the un-convolved ranges of templateflux
-        flux_lo = templateflux[:lo, :]
-        flux_hi = templateflux[hi:, :]
-        flux_mid = templateflux[lo:hi, :]
+        flux_lo = templateflux[:, :lo]
+        flux_hi = templateflux[:, hi:]
+        flux_mid = templateflux[:, lo:hi]
 
-        mid_len = flux_mid.shape[0]
+        mid_len = flux_mid.shape[1]
 
         fft_len = sc_fft.next_fast_len(mid_len + kernel_size - 1,
                                        real=True)
-        ft_flux_mid = sc_fft.rfft(flux_mid, n=fft_len, axis=0)
+        ft_flux_mid = sc_fft.rfft(flux_mid, n=fft_len)
 
-        return (np.vstack((flux_lo, flux_hi)), ft_flux_mid, fft_len)
+        return (np.hstack((flux_lo, flux_hi)), ft_flux_mid, fft_len)
 
 
     @staticmethod
-    def conv_pre_select(conv_pre, cols):
+    def conv_pre_select(conv_pre, rows):
         """
         Filter a set of preprocessing data for vdisp convolution
-        to use only the specified columns (templates)
+        to use only the specified rows (templates)
 
         Parameters
         ----------
         conv_pre: :class:`tuple` or None
             Preprocessing structure for vdisp convolution (may be None)
-        cols: :class:`int`
-            Which columns (templates) from the preprocessing
+        rows: :class:`int`
+            Which rows (templates) from the preprocessing
             data should we keep?
 
         Returns
         -------
-        Revised preprocessing data with only the selected columns
+        Revised preprocessing data with only the selected rows
 
         """
         if conv_pre is None:
             return None
         else:
             flux_lohi, ft_flux_mid, fft_len = conv_pre
-            return (flux_lohi[:, cols], ft_flux_mid[:, cols], fft_len)
+            return (flux_lohi[rows, :], ft_flux_mid[rows, :], fft_len)
 
 
     def convolve_vdisp_from_pre(self, flux_lohi, ft_flux_mid, flux_len, fft_len, vdisp):
@@ -322,7 +325,7 @@ def convolve_vdisp(self, templateflux, vdisp):
         ----------
         templateflux: :class:`np.ndarray`
            Either a 1D template array of fluxes, or a 2D array of size
-          [nfluxes x ntemplates] representing ntemplates fluxes
+          [ntemplates x nfluxes] representing ntemplates fluxes
         vdisp: :class:`float64`
            Velocity dispersion to convolve with fluxes
 
@@ -354,12 +357,12 @@ def convolve_vdisp(self, templateflux, vdisp):
                 output[:lo] = templateflux[:lo]
                 output[hi:] = templateflux[hi:]
             else:
-                ntemplates = templateflux.shape[1]
+                ntemplates = templateflux.shape[0]
                 for ii in range(ntemplates):
-                    output[lo:hi, ii] = self.convolve(
-                        templateflux[lo:hi, ii], kernel, mode='same')
-                output[:lo, :] = templateflux[:lo, :]
-                output[hi:, :] = templateflux[hi:, :]
+                    output[ii, lo:hi] = self.convolve(
+                        templateflux[ii, lo:hi], kernel, mode='same')
+                output[:, :lo] = templateflux[:, :lo]
+                output[:, hi:] = templateflux[:, hi:]
 
         return output