/cvx/sedumi/blkaux.c

  1/* ************************************************************
  2   MODULE sdmaux*.c  -- Several low-level subroutines for the
  3   mex-files in the Self-Dual-Minimization package.
  4
  5% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
  6% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
  7%
  8% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
  9%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
 10%   Supported by the Netherlands Organization for Scientific Research (NWO).
 11%
 12% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
 13%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
 14%
 15% Affiliations up to SeDuMi 1.02 (AUG1998):
 16%   CRL, McMaster University, Canada.
 17%   Supported by the Netherlands Organization for Scientific Research (NWO).
 18%
 19% This program is free software; you can redistribute it and/or modify
 20% it under the terms of the GNU General Public License as published by
 21% the Free Software Foundation; either version 2 of the License, or
 22% (at your option) any later version.
 23%
 24% This program is distributed in the hope that it will be useful,
 25% but WITHOUT ANY WARRANTY; without even the implied warranty of
 26% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 27% GNU General Public License for more details.
 28%
 29% You should have received a copy of the GNU General Public License
 30% along with this program; if not, write to the Free Software
 31% Foundation, Inc.,  51 Franklin Street, Fifth Floor, Boston, MA
 32% 02110-1301, USA
 33
 34   ************************************************************ */
 35
 36#include <string.h>
 37#include "blksdp.h"
 38
 39
 40/* ============================================================
 41   DOT-PRODUCT AND VECTOR ARRAY-OPS.
 42   ============================================================ */
 43
 44/* ************************************************************
 45   TIME-CRITICAL PROCEDURE -- realHadamard
 46   Computes  r = x .* y  using loop-unrolling.
 47   ************************************************************ */
 48void realHadamard(double * r, const double *x, const double *y, const mwIndex n)
 49{
 50/* TODO: This blas solution segfaults. OpenMP?
 51//     int one=1;
 52//     double zero=0.0;
 53//     #ifdef PC
 54//     dcopy(&n,y,&one,r,&one);
 55//     dtbmv('u','n','n',&n,&zero,x,&one,r,&one);
 56//     #endif
 57//     
 58//     #ifdef UNIX
 59//     dcopy_(&n,y,&one,r,&one);
 60//     dtbmv_('u','n','n',&n,&zero,x,&one,r,&one);
 61//     #endif
 62//     return;*/
 63    mwIndex i;
 64
 65 for(i=0; i+7< n; ){              /* LEVEL 8 */
 66   r[i] = x[i] * y[i]; i++;
 67   r[i] = x[i] * y[i]; i++;
 68   r[i] = x[i] * y[i]; i++;
 69   r[i] = x[i] * y[i]; i++;
 70   r[i] = x[i] * y[i]; i++;
 71   r[i] = x[i] * y[i]; i++;
 72   r[i] = x[i] * y[i]; i++;
 73   r[i] = x[i] * y[i]; i++;
 74 }
 75 if(i+3 < n){                         /* LEVEL 4 */
 76   r[i] = x[i] * y[i]; i++;
 77   r[i] = x[i] * y[i]; i++;
 78   r[i] = x[i] * y[i]; i++;
 79   r[i] = x[i] * y[i]; i++;
 80 }
 81/* ------------------------------------------------------------
 82   Now, i in {n-3, n-2, n-1, n}. Do the last n-i elements.
 83   ------------------------------------------------------------ */
 84 if(i+1 < n){                        /* LEVEL 2 */
 85   r[i] = x[i] * y[i]; i++;
 86   r[i] = x[i] * y[i]; i++;
 87 }
 88 if(i< n)                            /* LEVEL 1 */
 89   r[i] = x[i] * y[i];
 90}
 91
 92/* ************************************************************
 93   TIME-CRITICAL PROCEDURE -- realHadadiv
 94   Computes  r = x ./ y  using loop-unrolling.
 95   ************************************************************ */
 96void realHadadiv(double * r, const double *x, const double *y, const mwIndex n)
 97{
 98 mwIndex i;
 99
100 for(i=0; i+3< n; ){              /* LEVEL 4 */
101   r[i] = x[i] / y[i]; i++;
102   r[i] = x[i] / y[i]; i++;
103   r[i] = x[i] / y[i]; i++;
104   r[i] = x[i] / y[i]; i++;
105 }
106/* ------------------------------------------------------------
107   Now, i in {n-3, n-2, n-1, n}. Do the last n-i elements.
108   ------------------------------------------------------------ */
109 if(i+1 < n){                        /* LEVEL 2 */
110   r[i] = x[i] / y[i]; i++;
111   r[i] = x[i] / y[i]; i++;
112 }
113 if(i< n)                            /* LEVEL 1 */
114   r[i] = x[i] / y[i];
115}
116