/cvx/sedumi/blkaux.c

http://github.com/tcdoan/portfolioOptimization · C · 116 lines · 44 code · 8 blank · 64 comment · 7 complexity · 1dfbcba5e91523016fcaf29ef11505fd MD5 · raw file

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