Main Content

Generate Code to Optimize Portfolio by Using Black Litterman Approach

This example shows how to generate a MEX function and C source code from MATLAB® code that performs portfolio optimization using the Black Litterman approach.

Prerequisites

There are no prerequisites for this example.

About the hlblacklitterman Function

The hlblacklitterman.m function reads in financial information regarding a portfolio and performs portfolio optimization using the Black Litterman approach.

type hlblacklitterman
function [er, ps, w, pw, lambda, theta] = hlblacklitterman(delta, weq, sigma, tau, P, Q, Omega)%#codegen
% hlblacklitterman
%   This function performs the Black-Litterman blending of the prior
%   and the views into a new posterior estimate of the returns as
%   described in the paper by He and Litterman.
% Inputs
%   delta  - Risk tolerance from the equilibrium portfolio
%   weq    - Weights of the assets in the equilibrium portfolio
%   sigma  - Prior covariance matrix
%   tau    - Coefficiet of uncertainty in the prior estimate of the mean (pi)
%   P      - Pick matrix for the view(s)
%   Q      - Vector of view returns
%   Omega  - Matrix of variance of the views (diagonal)
% Outputs
%   Er     - Posterior estimate of the mean returns
%   w      - Unconstrained weights computed given the Posterior estimates
%            of the mean and covariance of returns.
%   lambda - A measure of the impact of each view on the posterior estimates.
%   theta  - A measure of the share of the prior and sample information in the
%            posterior precision.

% Reverse optimize and back out the equilibrium returns
% This is formula (12) page 6.
pi = weq * sigma * delta;
% We use tau * sigma many places so just compute it once
ts = tau * sigma;
% Compute posterior estimate of the mean
% This is a simplified version of formula (8) on page 4.
er = pi' + ts * P' * inv(P * ts * P' + Omega) * (Q - P * pi');
% We can also do it the long way to illustrate that d1 + d2 = I
d = inv(inv(ts) + P' * inv(Omega) * P);
d1 = d * inv(ts);
d2 = d * P' * inv(Omega) * P;
er2 = d1 * pi' + d2 * pinv(P) * Q;
% Compute posterior estimate of the uncertainty in the mean
% This is a simplified and combined version of formulas (9) and (15)
ps = ts - ts * P' * inv(P * ts * P' + Omega) * P * ts;
posteriorSigma = sigma + ps;
% Compute the share of the posterior precision from prior and views,
% then for each individual view so we can compare it with lambda
theta=zeros(1,2+size(P,1));
theta(1,1) = (trace(inv(ts) * ps) / size(ts,1));
theta(1,2) = (trace(P'*inv(Omega)*P* ps) / size(ts,1));
for i=1:size(P,1)
    theta(1,2+i) = (trace(P(i,:)'*inv(Omega(i,i))*P(i,:)* ps) / size(ts,1));
end
% Compute posterior weights based solely on changed covariance
w = (er' * inv(delta * posteriorSigma))';
% Compute posterior weights based on uncertainty in mean and covariance
pw = (pi * inv(delta * posteriorSigma))';
% Compute lambda value
% We solve for lambda from formula (17) page 7, rather than formula (18)
% just because it is less to type, and we've already computed w*.
lambda = pinv(P)' * (w'*(1+tau) - weq)';
end

% Black-Litterman example code for MatLab (hlblacklitterman.m)
% Copyright (c) Jay Walters, blacklitterman.org, 2008.
%
% Redistribution and use in source and binary forms, 
% with or without modification, are permitted provided 
% that the following conditions are met:
%
% Redistributions of source code must retain the above 
% copyright notice, this list of conditions and the following 
% disclaimer.
% 
% Redistributions in binary form must reproduce the above 
% copyright notice, this list of conditions and the following 
% disclaimer in the documentation and/or other materials 
% provided with the distribution.
%  
% Neither the name of blacklitterman.org nor the names of its
% contributors may be used to endorse or promote products 
% derived from this software without specific prior written
% permission.
%  
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND 
% CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, 
% INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 
% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, 
% WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 
% DAMAGE.
%
% This program uses the examples from the paper "The Intuition 
% Behind Black-Litterman Model  Portfolios", by He and Litterman,
% 1999.  You can find a copy of this  paper at the following url.
%     http:%papers.ssrn.com/sol3/papers.cfm?abstract_id=334304
%
% For more details on the Black-Litterman model you can also view
% "The BlackLitterman Model: A Detailed Exploration", by this author
% at the following url.
%     http:%www.blacklitterman.org/Black-Litterman.pdf
%

The %#codegen directive indicates that the MATLAB code is intended for code generation.

Generate the MEX Function for Testing

Generate a MEX function using the codegen command.

codegen hlblacklitterman -args {0, zeros(1, 7), zeros(7,7), 0, zeros(1, 7), 0, 0}
Code generation successful.

Before generating C code, you should first test the MEX function in MATLAB to ensure that it is functionally equivalent to the original MATLAB code and that no run-time errors occur. By default, codegen generates a MEX function named hlblacklitterman_mex in the current folder. This allows you to test the MATLAB code and MEX function and compare the results.

Run the MEX Function

Call the generated MEX function

testMex();
View 1
Country        P       mu      w*
Australia	     0	 4.328	 1.524
Canada   	     0	 7.576	 2.095
France   	 -29.5	 9.288	-3.948
Germany  	   100	 11.04	 35.41
Japan    	     0	 4.506	 11.05
UK       	 -70.5	 6.953	-9.462
USA      	     0	 8.069	 58.57
q        	     5
omega/tau	     0.0213
lambda   	     0.317
theta   	     0.0714
pr theta  	     0.929


View 1
Country        P       mu      w*
Australia	     0	 4.328	 1.524
Canada   	     0	 7.576	 2.095
France   	 -29.5	 9.288	-3.948
Germany  	   100	 11.04	 35.41
Japan    	     0	 4.506	 11.05
UK       	 -70.5	 6.953	-9.462
USA      	     0	 8.069	 58.57
q        	     5
omega/tau	     0.0213
lambda   	     0.317
theta   	     0.0714
pr theta  	     0.929

Execution Time - MATLAB function: 0.034639 seconds
Execution Time - MEX function   : 0.013608 seconds

Generate C Code

cfg = coder.config('lib');
codegen -config cfg hlblacklitterman  -args {0, zeros(1, 7), zeros(7,7), 0, zeros(1, 7), 0, 0}
Code generation successful.

Using codegen with the specified -config cfg option produces a standalone C library.

Inspect the Generated Code

By default, the code generated for the library is in the folder codegen/lib/hbblacklitterman/.

The files are:

dir codegen/lib/hlblacklitterman/
.                              ..                             _clang-format                  buildInfo.mat                  codeInfo.mat                   codedescriptor.dmr             compileInfo.mat                examples                       hlblacklitterman.a             hlblacklitterman.c             hlblacklitterman.h             hlblacklitterman.o             hlblacklitterman_data.h        hlblacklitterman_initialize.c  hlblacklitterman_initialize.h  hlblacklitterman_initialize.o  hlblacklitterman_rtw.mk        hlblacklitterman_terminate.c   hlblacklitterman_terminate.h   hlblacklitterman_terminate.o   hlblacklitterman_types.h       interface                      inv.c                          inv.h                          inv.o                          rtGetInf.c                     rtGetInf.h                     rtGetInf.o                     rtGetNaN.c                     rtGetNaN.h                     rtGetNaN.o                     rt_nonfinite.c                 rt_nonfinite.h                 rt_nonfinite.o                 rtw_proj.tmw                   rtwtypes.h                     

Inspect the C Code for the hlblacklitterman.c Function

type codegen/lib/hlblacklitterman/hlblacklitterman.c
/*
 * File: hlblacklitterman.c
 *
 * MATLAB Coder version            : 24.1
 * C/C++ source code generated on  : 12-Feb-2024 21:03:37
 */

/* Include Files */
#include "hlblacklitterman.h"
#include "inv.h"
#include "rt_nonfinite.h"
#include "rt_nonfinite.h"
#include <emmintrin.h>
#include <math.h>

/* Function Definitions */
/*
 * hlblacklitterman
 *    This function performs the Black-Litterman blending of the prior
 *    and the views into a new posterior estimate of the returns as
 *    described in the paper by He and Litterman.
 *  Inputs
 *    delta  - Risk tolerance from the equilibrium portfolio
 *    weq    - Weights of the assets in the equilibrium portfolio
 *    sigma  - Prior covariance matrix
 *    tau    - Coefficiet of uncertainty in the prior estimate of the mean (pi)
 *    P      - Pick matrix for the view(s)
 *    Q      - Vector of view returns
 *    Omega  - Matrix of variance of the views (diagonal)
 *  Outputs
 *    Er     - Posterior estimate of the mean returns
 *    w      - Unconstrained weights computed given the Posterior estimates
 *             of the mean and covariance of returns.
 *    lambda - A measure of the impact of each view on the posterior estimates.
 *    theta  - A measure of the share of the prior and sample information in the
 *             posterior precision.
 *
 * Arguments    : double delta
 *                const double weq[7]
 *                const double sigma[49]
 *                double tau
 *                const double P[7]
 *                double Q
 *                double Omega
 *                double er[7]
 *                double ps[49]
 *                double w[7]
 *                double pw[7]
 *                double *lambda
 *                double theta[3]
 * Return Type  : void
 */
void hlblacklitterman(double delta, const double weq[7], const double sigma[49],
                      double tau, const double P[7], double Q, double Omega,
                      double er[7], double ps[49], double w[7], double pw[7],
                      double *lambda, double theta[3])
{
  __m128d r;
  __m128d r1;
  double b_er_tmp[49];
  double dv[49];
  double posteriorSigma[49];
  double ts[49];
  double er_tmp[7];
  double pi[7];
  double y_tmp[7];
  double absxk;
  double anrm;
  double cfromc;
  double ctoc;
  double nrm;
  double scale;
  int br;
  int cr;
  int exponent;
  int ib;
  int ic;
  int scalarLB;
  boolean_T doscale;
  /*  Reverse optimize and back out the equilibrium returns */
  /*  This is formula (12) page 6. */
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    scale = 0.0;
    for (ib = 0; ib < 7; ib++) {
      scale += weq[ib] * sigma[ib + 7 * scalarLB];
    }
    pi[scalarLB] = scale * delta;
  }
  /*  We use tau * sigma many places so just compute it once */
  for (scalarLB = 0; scalarLB <= 46; scalarLB += 2) {
    _mm_storeu_pd(&ts[scalarLB],
                  _mm_mul_pd(_mm_set1_pd(tau), _mm_loadu_pd(&sigma[scalarLB])));
  }
  ts[48] = tau * sigma[48];
  /*  Compute posterior estimate of the mean */
  /*  This is a simplified version of formula (8) on page 4. */
  anrm = 0.0;
  cfromc = 0.0;
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    scale = 0.0;
    absxk = 0.0;
    for (ib = 0; ib < 7; ib++) {
      nrm = P[ib];
      scale += ts[scalarLB + 7 * ib] * nrm;
      absxk += nrm * ts[ib + 7 * scalarLB];
    }
    y_tmp[scalarLB] = absxk;
    er_tmp[scalarLB] = scale;
    scale = P[scalarLB];
    anrm += absxk * scale;
    cfromc += scale * pi[scalarLB];
  }
  ctoc = 1.0 / (anrm + Omega);
  absxk = Q - cfromc;
  /*  We can also do it the long way to illustrate that d1 + d2 = I */
  anrm = 1.0 / Omega;
  /*  Compute posterior estimate of the uncertainty in the mean */
  /*  This is a simplified and combined version of formulas (9) and (15) */
  nrm = 0.0;
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    er[scalarLB] = pi[scalarLB] + er_tmp[scalarLB] * ctoc * absxk;
    nrm += y_tmp[scalarLB] * P[scalarLB];
  }
  ctoc = 1.0 / (nrm + Omega);
  r = _mm_set1_pd(ctoc);
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    __m128d r2;
    r1 = _mm_loadu_pd(&er_tmp[0]);
    scale = P[scalarLB];
    r2 = _mm_set1_pd(scale);
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB], _mm_mul_pd(_mm_mul_pd(r1, r), r2));
    r1 = _mm_loadu_pd(&er_tmp[2]);
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB + 2],
                  _mm_mul_pd(_mm_mul_pd(r1, r), r2));
    r1 = _mm_loadu_pd(&er_tmp[4]);
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB + 4],
                  _mm_mul_pd(_mm_mul_pd(r1, r), r2));
    b_er_tmp[7 * scalarLB + 6] = er_tmp[6] * ctoc * scale;
  }
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    for (ib = 0; ib < 7; ib++) {
      scale = 0.0;
      for (br = 0; br < 7; br++) {
        scale += b_er_tmp[scalarLB + 7 * br] * ts[br + 7 * ib];
      }
      br = scalarLB + 7 * ib;
      ps[br] = ts[br] - scale;
    }
  }
  for (scalarLB = 0; scalarLB <= 46; scalarLB += 2) {
    r = _mm_loadu_pd(&ps[scalarLB]);
    _mm_storeu_pd(&posteriorSigma[scalarLB],
                  _mm_add_pd(_mm_loadu_pd(&sigma[scalarLB]), r));
  }
  posteriorSigma[48] = sigma[48] + ps[48];
  /*  Compute the share of the posterior precision from prior and views, */
  /*  then for each individual view so we can compare it with lambda */
  inv(ts, dv);
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    for (ib = 0; ib < 7; ib++) {
      scale = 0.0;
      for (br = 0; br < 7; br++) {
        scale += dv[scalarLB + 7 * br] * ps[br + 7 * ib];
      }
      ts[scalarLB + 7 * ib] = scale;
    }
  }
  cfromc = 0.0;
  for (br = 0; br < 7; br++) {
    cfromc += ts[br + 7 * br];
  }
  theta[0] = cfromc / 7.0;
  r = _mm_set1_pd(anrm);
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    r1 = _mm_set1_pd(P[scalarLB]);
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB],
                  _mm_mul_pd(_mm_mul_pd(_mm_loadu_pd(&P[0]), r), r1));
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB + 2],
                  _mm_mul_pd(_mm_mul_pd(_mm_loadu_pd(&P[2]), r), r1));
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB + 4],
                  _mm_mul_pd(_mm_mul_pd(_mm_loadu_pd(&P[4]), r), r1));
    b_er_tmp[7 * scalarLB + 6] = P[6] * anrm * P[scalarLB];
  }
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    for (ib = 0; ib < 7; ib++) {
      scale = 0.0;
      for (br = 0; br < 7; br++) {
        scale += b_er_tmp[scalarLB + 7 * br] * ps[br + 7 * ib];
      }
      ts[scalarLB + 7 * ib] = scale;
    }
  }
  cfromc = 0.0;
  for (br = 0; br < 7; br++) {
    cfromc += ts[br + 7 * br];
  }
  theta[1] = cfromc / 7.0;
  r = _mm_set1_pd(anrm);
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    r1 = _mm_set1_pd(P[scalarLB]);
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB],
                  _mm_mul_pd(_mm_mul_pd(_mm_loadu_pd(&P[0]), r), r1));
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB + 2],
                  _mm_mul_pd(_mm_mul_pd(_mm_loadu_pd(&P[2]), r), r1));
    _mm_storeu_pd(&b_er_tmp[7 * scalarLB + 4],
                  _mm_mul_pd(_mm_mul_pd(_mm_loadu_pd(&P[4]), r), r1));
    b_er_tmp[7 * scalarLB + 6] = P[6] * anrm * P[scalarLB];
  }
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    for (ib = 0; ib < 7; ib++) {
      scale = 0.0;
      for (br = 0; br < 7; br++) {
        scale += b_er_tmp[scalarLB + 7 * br] * ps[br + 7 * ib];
      }
      ts[scalarLB + 7 * ib] = scale;
    }
  }
  cfromc = 0.0;
  for (br = 0; br < 7; br++) {
    cfromc += ts[br + 7 * br];
  }
  theta[2] = cfromc / 7.0;
  /*  Compute posterior weights based solely on changed covariance */
  for (scalarLB = 0; scalarLB <= 46; scalarLB += 2) {
    r = _mm_loadu_pd(&posteriorSigma[scalarLB]);
    _mm_storeu_pd(&b_er_tmp[scalarLB], _mm_mul_pd(_mm_set1_pd(delta), r));
  }
  b_er_tmp[48] = delta * posteriorSigma[48];
  inv(b_er_tmp, dv);
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    scale = 0.0;
    for (ib = 0; ib < 7; ib++) {
      scale += er[ib] * dv[ib + 7 * scalarLB];
    }
    w[scalarLB] = scale;
  }
  /*  Compute posterior weights based on uncertainty in mean and covariance */
  for (scalarLB = 0; scalarLB <= 46; scalarLB += 2) {
    r = _mm_loadu_pd(&posteriorSigma[scalarLB]);
    _mm_storeu_pd(&posteriorSigma[scalarLB], _mm_mul_pd(_mm_set1_pd(delta), r));
  }
  posteriorSigma[48] *= delta;
  inv(posteriorSigma, dv);
  /*  Compute lambda value */
  /*  We solve for lambda from formula (17) page 7, rather than formula (18) */
  /*  just because it is less to type, and we've already computed w*. */
  for (scalarLB = 0; scalarLB < 7; scalarLB++) {
    scale = 0.0;
    for (ib = 0; ib < 7; ib++) {
      scale += pi[ib] * dv[ib + 7 * scalarLB];
    }
    pw[scalarLB] = scale;
    er_tmp[scalarLB] = P[scalarLB];
  }
  doscale = true;
  for (br = 0; br < 7; br++) {
    pi[br] = 0.0;
    if (doscale) {
      scale = P[br];
      if (rtIsInf(scale) || rtIsNaN(scale)) {
        doscale = false;
      }
    } else {
      doscale = false;
    }
  }
  if (!doscale) {
    for (scalarLB = 0; scalarLB < 7; scalarLB++) {
      pi[scalarLB] = rtNaN;
    }
  } else {
    double cscale;
    boolean_T exitg1;
    boolean_T guard1;
    boolean_T notdone;
    doscale = false;
    anrm = 0.0;
    br = 0;
    exitg1 = false;
    while ((!exitg1) && (br < 7)) {
      absxk = fabs(P[br]);
      if (rtIsNaN(absxk)) {
        anrm = rtNaN;
        exitg1 = true;
      } else {
        if (absxk > anrm) {
          anrm = absxk;
        }
        br++;
      }
    }
    cscale = anrm;
    guard1 = false;
    if ((anrm > 0.0) && (anrm < 6.7178761075670888E-139)) {
      doscale = true;
      cscale = 6.7178761075670888E-139;
      guard1 = true;
    } else if (anrm > 1.4885657073574029E+138) {
      doscale = true;
      cscale = 1.4885657073574029E+138;
      guard1 = true;
    }
    if (guard1) {
      cfromc = anrm;
      ctoc = cscale;
      notdone = true;
      while (notdone) {
        absxk = cfromc * 2.0041683600089728E-292;
        nrm = ctoc / 4.9896007738368E+291;
        if ((absxk > ctoc) && (ctoc != 0.0)) {
          scale = 2.0041683600089728E-292;
          cfromc = absxk;
        } else if (nrm > cfromc) {
          scale = 4.9896007738368E+291;
          ctoc = nrm;
        } else {
          scale = ctoc / cfromc;
          notdone = false;
        }
        r = _mm_loadu_pd(&er_tmp[0]);
        r1 = _mm_set1_pd(scale);
        _mm_storeu_pd(&er_tmp[0], _mm_mul_pd(r, r1));
        r = _mm_loadu_pd(&er_tmp[2]);
        _mm_storeu_pd(&er_tmp[2], _mm_mul_pd(r, r1));
        r = _mm_loadu_pd(&er_tmp[4]);
        _mm_storeu_pd(&er_tmp[4], _mm_mul_pd(r, r1));
        er_tmp[6] *= scale;
      }
    }
    nrm = 0.0;
    scale = 3.3121686421112381E-170;
    for (br = 0; br < 7; br++) {
      absxk = fabs(er_tmp[br]);
      if (absxk > scale) {
        cfromc = scale / absxk;
        nrm = nrm * cfromc * cfromc + 1.0;
        scale = absxk;
      } else {
        cfromc = absxk / scale;
        nrm += cfromc * cfromc;
      }
    }
    nrm = scale * sqrt(nrm);
    if (nrm > 0.0) {
      if (er_tmp[0] < 0.0) {
        cfromc = -nrm;
      } else {
        cfromc = nrm;
      }
      if (fabs(cfromc) >= 1.0020841800044864E-292) {
        absxk = 1.0 / cfromc;
        r = _mm_loadu_pd(&er_tmp[0]);
        r1 = _mm_set1_pd(absxk);
        _mm_storeu_pd(&er_tmp[0], _mm_mul_pd(r1, r));
        r = _mm_loadu_pd(&er_tmp[2]);
        _mm_storeu_pd(&er_tmp[2], _mm_mul_pd(r1, r));
        r = _mm_loadu_pd(&er_tmp[4]);
        _mm_storeu_pd(&er_tmp[4], _mm_mul_pd(r1, r));
        er_tmp[6] *= absxk;
      } else {
        r = _mm_loadu_pd(&er_tmp[0]);
        r1 = _mm_set1_pd(cfromc);
        _mm_storeu_pd(&er_tmp[0], _mm_div_pd(r, r1));
        r = _mm_loadu_pd(&er_tmp[2]);
        _mm_storeu_pd(&er_tmp[2], _mm_div_pd(r, r1));
        r = _mm_loadu_pd(&er_tmp[4]);
        _mm_storeu_pd(&er_tmp[4], _mm_div_pd(r, r1));
        er_tmp[6] /= cfromc;
      }
      er_tmp[0]++;
      cfromc = -cfromc;
    } else {
      cfromc = 0.0;
    }
    for (br = 0; br < 7; br++) {
      y_tmp[br] = er_tmp[br];
    }
    if (cfromc != 0.0) {
      r = _mm_loadu_pd(&y_tmp[0]);
      r1 = _mm_set1_pd(-1.0);
      _mm_storeu_pd(&y_tmp[0], _mm_mul_pd(r, r1));
      r = _mm_loadu_pd(&y_tmp[2]);
      _mm_storeu_pd(&y_tmp[2], _mm_mul_pd(r, r1));
      r = _mm_loadu_pd(&y_tmp[4]);
      _mm_storeu_pd(&y_tmp[4], _mm_mul_pd(r, r1));
      y_tmp[6] = -y_tmp[6];
      y_tmp[0]++;
      nrm = fabs(cfromc);
      absxk = cfromc / nrm;
      cfromc = nrm;
      r = _mm_loadu_pd(&y_tmp[0]);
      r1 = _mm_set1_pd(absxk);
      _mm_storeu_pd(&y_tmp[0], _mm_mul_pd(r1, r));
      r = _mm_loadu_pd(&y_tmp[2]);
      _mm_storeu_pd(&y_tmp[2], _mm_mul_pd(r1, r));
      r = _mm_loadu_pd(&y_tmp[4]);
      _mm_storeu_pd(&y_tmp[4], _mm_mul_pd(r1, r));
      y_tmp[6] *= absxk;
    } else {
      for (br = 0; br < 7; br++) {
        y_tmp[br] = 0.0;
      }
      y_tmp[0] = 1.0;
    }
    if (doscale) {
      notdone = true;
      while (notdone) {
        absxk = cscale * 2.0041683600089728E-292;
        nrm = anrm / 4.9896007738368E+291;
        if ((absxk > anrm) && (anrm != 0.0)) {
          scale = 2.0041683600089728E-292;
          cscale = absxk;
        } else if (nrm > cscale) {
          scale = 4.9896007738368E+291;
          anrm = nrm;
        } else {
          scale = anrm / cscale;
          notdone = false;
        }
        cfromc *= scale;
      }
    }
    doscale = (rtIsInf(cfromc) || rtIsNaN(cfromc));
    if (doscale) {
      absxk = rtNaN;
    } else if (cfromc < 4.4501477170144028E-308) {
      absxk = 4.94065645841247E-324;
    } else {
      frexp(cfromc, &exponent);
      absxk = ldexp(1.0, exponent - 53);
    }
    absxk *= 7.0;
    if (doscale) {
      absxk = 1.7976931348623157E+308;
    }
    if (cfromc > absxk) {
      absxk = 1.0 / cfromc;
      for (cr = 0; cr < 7; cr++) {
        scalarLB = cr + 1;
        for (ic = scalarLB; ic <= scalarLB; ic++) {
          pi[ic - 1] = 0.0;
        }
      }
      br = 0;
      for (cr = 0; cr < 7; cr++) {
        br++;
        for (ib = br; ib <= br; ib += 7) {
          scalarLB = cr + 1;
          exponent = cr - 1;
          for (ic = scalarLB; ic <= exponent; ic += 2) {
            r = _mm_loadu_pd(&pi[ic - 1]);
            _mm_storeu_pd(&pi[ic - 1],
                          _mm_add_pd(r, _mm_set1_pd(y_tmp[ib - 1] * absxk)));
          }
          for (ic = scalarLB; ic <= scalarLB; ic++) {
            pi[ic - 1] += y_tmp[ib - 1] * absxk;
          }
        }
      }
    }
  }
  *lambda = 0.0;
  for (br = 0; br < 7; br++) {
    *lambda += pi[br] * (w[br] * (tau + 1.0) - weq[br]);
  }
}

/*
 * File trailer for hlblacklitterman.c
 *
 * [EOF]
 */