/*
 * mvnrnd_.c
 *
 * Code generation for function 'mvnrnd_'
 *
 * C source code generated on: Wed Aug 26 14:59:33 2015
 *
 */

/* Include files */
#include "rt_nonfinite.h"
#include "Select_Ground_Motions.h"
#include "mvnrnd_.h"
#include "Select_Ground_Motions_emxutil.h"
#include "inv.h"
#include "sqrt.h"
#include "cholcov_.h"
#include "Select_Ground_Motions_mexutil.h"
#include "Select_Ground_Motions_data.h"

/* Variable Definitions */
static emlrtRSInfo cf_emlrtRSI = { 100, "mvnrnd_",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/mvnrnd_.m"
};

static emlrtRSInfo df_emlrtRSI = { 104, "mvnrnd_",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/mvnrnd_.m"
};

static emlrtRSInfo ef_emlrtRSI = { 110, "mvnrnd_",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/mvnrnd_.m"
};

static emlrtRSInfo ff_emlrtRSI = { 114, "mvnrnd_",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/mvnrnd_.m"
};

static emlrtRSInfo gf_emlrtRSI = { 88, "randn",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/randfun/randn.m" };

static emlrtRSInfo hf_emlrtRSI = { 21, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtRSInfo if_emlrtRSI = { 23, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtRSInfo jf_emlrtRSI = { 72, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtRSInfo kf_emlrtRSI = { 78, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtRSInfo lf_emlrtRSI = { 84, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtMCInfo db_emlrtMCI = { 88, 9, "randn",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/randfun/randn.m" };

static emlrtMCInfo eb_emlrtMCI = { 22, 5, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtMCInfo fb_emlrtMCI = { 21, 15, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtMCInfo gb_emlrtMCI = { 24, 5, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtMCInfo hb_emlrtMCI = { 23, 15, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtRTEInfo fb_emlrtRTEI = { 1, 18, "mvnrnd_",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/mvnrnd_.m"
};

static emlrtRTEInfo gb_emlrtRTEI = { 60, 1, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtRTEInfo hb_emlrtRTEI = { 63, 1, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtRTEInfo ib_emlrtRTEI = { 100, 9, "mvnrnd_",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/mvnrnd_.m"
};

static emlrtRTEInfo jb_emlrtRTEI = { 82, 5, "bsxfun",
  "/usr/local/MATLAB/R2013a/toolbox/eml/lib/matlab/elmat/bsxfun.m" };

static emlrtECInfo cb_emlrtECI = { 2, 104, 13, "mvnrnd_",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/mvnrnd_.m"
};

static emlrtECInfo db_emlrtECI = { 2, 114, 13, "mvnrnd_",
  "/home/sebenik/Projects/Ground_Motion_Selection/GMS_con_cms_1gmc_dynamicInput_501_jure/mvnrnd_.m"
};

/* Function Declarations */
static void c_emlrt_marshallIn(const mxArray *b_randn, const char_T *identifier,
  emxArray_real_T *y);
static void d_emlrt_marshallIn(const mxArray *u, const emlrtMsgIdentifier
  *parentId, emxArray_real_T *y);
static void l_emlrt_marshallIn(const mxArray *src, const emlrtMsgIdentifier
  *msgId, emxArray_real_T *ret);
static const mxArray *randn(const mxArray *b, const mxArray *c, emlrtMCInfo
  *location);

/* Function Definitions */
static void c_emlrt_marshallIn(const mxArray *b_randn, const char_T *identifier,
  emxArray_real_T *y)
{
  emlrtMsgIdentifier thisId;
  thisId.fIdentifier = identifier;
  thisId.fParent = NULL;
  d_emlrt_marshallIn(emlrtAlias(b_randn), &thisId, y);
  emlrtDestroyArray(&b_randn);
}

static void d_emlrt_marshallIn(const mxArray *u, const emlrtMsgIdentifier
  *parentId, emxArray_real_T *y)
{
  l_emlrt_marshallIn(emlrtAlias(u), parentId, y);
  emlrtDestroyArray(&u);
}

static void l_emlrt_marshallIn(const mxArray *src, const emlrtMsgIdentifier
  *msgId, emxArray_real_T *ret)
{
  int32_T iv135[2];
  boolean_T bv1[2];
  int32_T i25;
  static const boolean_T bv2[2] = { FALSE, TRUE };

  int32_T iv136[2];
  for (i25 = 0; i25 < 2; i25++) {
    iv135[i25] = 1 + -2 * i25;
    bv1[i25] = bv2[i25];
  }

  emlrtCheckVsBuiltInR2012b(emlrtRootTLSGlobal, msgId, src, "double", FALSE, 2U,
    iv135, bv1, iv136);
  i25 = ret->size[0] * ret->size[1];
  ret->size[0] = iv136[0];
  ret->size[1] = iv136[1];
  emxEnsureCapacity((emxArray__common *)ret, i25, (int32_T)sizeof(real_T),
                    (emlrtRTEInfo *)NULL);
  emlrtImportArrayR2011b(src, ret->data, 8, FALSE);
  emlrtDestroyArray(&src);
}

static const mxArray *randn(const mxArray *b, const mxArray *c, emlrtMCInfo
  *location)
{
  const mxArray *pArrays[2];
  const mxArray *m55;
  pArrays[0] = b;
  pArrays[1] = c;
  return emlrtCallMATLABR2012b(emlrtRootTLSGlobal, 1, &m55, 2, pArrays, "randn",
    TRUE, location);
}

void mvnrnd_(const emxArray_real_T *mu, const emxArray_real_T *sigma,
             emxArray_real_T *r)
{
  uint32_T sz[2];
  int32_T i6;
  boolean_T sigmaIsDiag;
  emxArray_real_T *t;
  emxArray_real_T *b_r;
  int32_T loop_ub;
  const mxArray *y;
  static const int32_T iv32[2] = { 1, 44 };

  const mxArray *m8;
  char_T cv34[44];
  int32_T i;
  static const char_T cv35[44] = { 'C', 'o', 'd', 'e', 'r', ':', 'M', 'A', 'T',
    'L', 'A', 'B', ':', 'b', 's', 'x', 'f', 'u', 'n', '_', 'a', 'r', 'r', 'a',
    'y', 'D', 'i', 'm', 'e', 'n', 's', 'i', 'o', 'n', 's', 'M', 'u', 's', 't',
    'M', 'a', 't', 'c', 'h' };

  const mxArray *b_y;
  static const int32_T iv33[2] = { 1, 37 };

  char_T cv36[37];
  static const char_T cv37[37] = { 'C', 'o', 'd', 'e', 'r', ':', 't', 'o', 'o',
    'l', 'b', 'o', 'x', ':', 'b', 's', 'x', 'f', 'u', 'n', '_', 'd', 'y', 'n',
    'a', 'm', 'i', 'c', 'E', 'x', 'p', 'a', 'n', 's', 'i', 'o', 'n' };

  emxArray_real_T *r22;
  emxArray_real_T *bv;
  int32_T asub;
  int32_T bsub;
  int32_T ak;
  int32_T bk;
  int32_T c_r;
  int32_T ck;
  emxArray_real_T *cv;
  int32_T exitg1;
  real_T err;
  int32_T exitg2;
  int32_T d_r[2];
  int32_T b_mu[2];
  const mxArray *c_y;
  static const int32_T iv34[2] = { 1, 45 };

  char_T cv38[45];
  static const char_T cv39[45] = { 'C', 'o', 'd', 'e', 'r', ':', 't', 'o', 'o',
    'l', 'b', 'o', 'x', ':', 'm', 't', 'i', 'm', 'e', 's', '_', 'n', 'o', 'D',
    'y', 'n', 'a', 'm', 'i', 'c', 'S', 'c', 'a', 'l', 'a', 'r', 'E', 'x', 'p',
    'a', 'n', 's', 'i', 'o', 'n' };

  const mxArray *d_y;
  static const int32_T iv35[2] = { 1, 21 };

  char_T cv40[21];
  static const char_T cv41[21] = { 'C', 'o', 'd', 'e', 'r', ':', 'M', 'A', 'T',
    'L', 'A', 'B', ':', 'i', 'n', 'n', 'e', 'r', 'd', 'i', 'm' };

  emxArray_real_T *r23;
  real_T beta1;
  char_T TRANSB;
  char_T TRANSA;
  ptrdiff_t m_t;
  ptrdiff_t n_t;
  ptrdiff_t k_t;
  ptrdiff_t lda_t;
  ptrdiff_t ldb_t;
  ptrdiff_t ldc_t;
  double * alpha1_t;
  double * Aia0_t;
  double * Bib0_t;
  double * beta1_t;
  double * Cic0_t;
  emlrtHeapReferenceStackEnterFcnR2012b(emlrtRootTLSGlobal);

  /* MVNRND Random vectors from the multivariate normal distribution. */
  /*    R = MVNRND(MU,SIGMA) returns an N-by-D matrix R of random vectors */
  /*    chosen from the multivariate normal distribution with mean vector MU, */
  /*    and covariance matrix SIGMA.  MU is an N-by-D matrix, and MVNRND */
  /*    generates each row of R using the corresponding row of MU.  SIGMA is a */
  /*    D-by-D symmetric positive semi-definite matrix, or a D-by-D-by-N array. */
  /*    If SIGMA is an array, MVNRND generates each row of R using the */
  /*    corresponding page of SIGMA, i.e., MVNRND computes R(I,:) using MU(I,:) */
  /*    and SIGMA(:,:,I).  If the covariance matrix is diagonal, containing */
  /*    variances along the diagonal and zero covariances off the diagonal, */
  /*    SIGMA may also be specified as a 1-by-D matrix or a 1-by-D-by-N array, */
  /*    containing just the diagonal. If MU is a 1-by-D vector, MVNRND */
  /*    replicates it to match the trailing dimension of SIGMA. */
  /*  */
  /*    R = MVNRND(MU,SIGMA,N) returns a N-by-D matrix R of random vectors */
  /*    chosen from the multivariate normal distribution with 1-by-D mean */
  /*    vector MU, and D-by-D covariance matrix SIGMA. */
  /*  */
  /*    Example: */
  /*  */
  /*       mu = [1 -1]; Sigma = [.9 .4; .4 .3]; */
  /*       r = mvnrnd(mu, Sigma, 500); */
  /*       plot(r(:,1),r(:,2),'.'); */
  /*  */
  /*    See also MVTRND, MVNPDF, MVNCDF, NORMRND. */
  /*    R = MVNRND(MU,SIGMA,N,T) supplies the Cholesky factor T of */
  /*    SIGMA, so that SIGMA(:,:,J) == T(:,:,J)'*T(:,:,J) if SIGMA is a 3D array or SIGMA */
  /*    == T'*T if SIGMA is a matrix.  No error checking is done on T. */
  /*  */
  /*    [R,T] = MVNRND(...) returns the Cholesky factor T, so it can be */
  /*    re-used to make later calls more efficient, although there are greater */
  /*    efficiency gains when SIGMA can be specified as a diagonal instead. */
  /*    Copyright 1993-2011 The MathWorks, Inc. */
  for (i6 = 0; i6 < 2; i6++) {
    sz[i6] = (uint32_T)sigma->size[i6];
  }

  if ((sz[0] == 1U) && (sz[1] > 1U)) {
    /*  Just the diagonal of Sigma has been passed in. */
    sigmaIsDiag = TRUE;
  } else {
    sigmaIsDiag = FALSE;
  }

  /*  Special case: if mu is a column vector, then use sigma to try */
  /*  to interpret it as a row vector. */
  /* if d == 1 && sz(1) == n */
  /*     mu = mu.'; */
  /*     [n,d] = size(mu); */
  /* end */
  /*  Get size of data. */
  /*  cases not supplied */
  /*  Single covariance matrix */
  /*  Make sure sigma is the right size */
  b_emxInit_real_T(&t, 2, &ib_emlrtRTEI, TRUE);
  b_emxInit_real_T(&b_r, 2, &fb_emlrtRTEI, TRUE);
  if (sigmaIsDiag) {
    /*  Just the diagonal of sigma has been specified. */
    /* if not(all(sigma>=0)) */
    /* error(message('stats:mvnrnd:BadDiagSigma')); */
    /* end */
    emlrtPushRtStackR2012b(&cf_emlrtRSI, emlrtRootTLSGlobal);
    i6 = t->size[0] * t->size[1];
    t->size[0] = sigma->size[0];
    t->size[1] = sigma->size[1];
    emxEnsureCapacity((emxArray__common *)t, i6, (int32_T)sizeof(real_T),
                      &fb_emlrtRTEI);
    loop_ub = sigma->size[0] * sigma->size[1];
    for (i6 = 0; i6 < loop_ub; i6++) {
      t->data[i6] = sigma->data[i6];
    }

    d_sqrt(t);
    emlrtPopRtStackR2012b(&cf_emlrtRSI, emlrtRootTLSGlobal);
    emlrtPushRtStackR2012b(&df_emlrtRSI, emlrtRootTLSGlobal);
    emlrtPushRtStackR2012b(&gf_emlrtRSI, emlrtRootTLSGlobal);
    c_emlrt_marshallIn(randn(emlrt_marshallOut(1.0), emlrt_marshallOut(mu->size
      [1]), &db_emlrtMCI), "randn", b_r);
    emlrtPopRtStackR2012b(&gf_emlrtRSI, emlrtRootTLSGlobal);
    if ((b_r->size[1] != 1) && (t->size[1] != 1) && (b_r->size[1] != t->size[1]))
    {
      sigmaIsDiag = FALSE;
    } else {
      sigmaIsDiag = TRUE;
    }

    if (sigmaIsDiag) {
    } else {
      emlrtPushRtStackR2012b(&hf_emlrtRSI, emlrtRootTLSGlobal);
      y = NULL;
      m8 = mxCreateCharArray(2, iv32);
      for (i = 0; i < 44; i++) {
        cv34[i] = cv35[i];
      }

      emlrtInitCharArrayR2013a(emlrtRootTLSGlobal, 44, m8, cv34);
      emlrtAssign(&y, m8);
      error(message(y, &eb_emlrtMCI), &fb_emlrtMCI);
      emlrtPopRtStackR2012b(&hf_emlrtRSI, emlrtRootTLSGlobal);
    }

    if (b_r->size[1] != t->size[1]) {
      sigmaIsDiag = FALSE;
    } else {
      sigmaIsDiag = TRUE;
    }

    if (sigmaIsDiag) {
    } else {
      emlrtPushRtStackR2012b(&if_emlrtRSI, emlrtRootTLSGlobal);
      b_y = NULL;
      m8 = mxCreateCharArray(2, iv33);
      for (i = 0; i < 37; i++) {
        cv36[i] = cv37[i];
      }

      emlrtInitCharArrayR2013a(emlrtRootTLSGlobal, 37, m8, cv36);
      emlrtAssign(&b_y, m8);
      error(message(b_y, &gb_emlrtMCI), &hb_emlrtMCI);
      emlrtPopRtStackR2012b(&if_emlrtRSI, emlrtRootTLSGlobal);
    }

    b_emxInit_real_T(&r22, 2, &fb_emlrtRTEI, TRUE);
    emxInit_real_T(&bv, 1, &hb_emlrtRTEI, TRUE);
    i6 = r22->size[0] * r22->size[1];
    r22->size[0] = t->size[0];
    r22->size[1] = b_r->size[1];
    emxEnsureCapacity((emxArray__common *)r22, i6, (int32_T)sizeof(real_T),
                      &gb_emlrtRTEI);
    i = t->size[0];
    i6 = bv->size[0];
    bv->size[0] = i;
    emxEnsureCapacity((emxArray__common *)bv, i6, (int32_T)sizeof(real_T),
                      &hb_emlrtRTEI);
    asub = 1;
    bsub = 1;
    ak = 0;
    bk = -1;
    loop_ub = t->size[0];
    i = b_r->size[1];
    c_r = t->size[0];
    c_r = loop_ub * i - c_r;
    emlrtPushRtStackR2012b(&jf_emlrtRSI, emlrtRootTLSGlobal);
    loop_ub = t->size[0];
    if ((loop_ub == 0) || (0 > c_r)) {
      sigmaIsDiag = FALSE;
    } else {
      loop_ub = t->size[0];
      sigmaIsDiag = (c_r > MAX_int32_T - loop_ub);
    }

    if (sigmaIsDiag) {
      emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
      check_forloop_overflow_error();
      emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
    }

    emlrtPopRtStackR2012b(&jf_emlrtRSI, emlrtRootTLSGlobal);
    ck = 0;
    emxInit_real_T(&cv, 1, &jb_emlrtRTEI, TRUE);
    do {
      exitg1 = 0;
      loop_ub = t->size[0];
      if ((loop_ub > 0) && (ck <= c_r)) {
        emlrtPushRtStackR2012b(&kf_emlrtRSI, emlrtRootTLSGlobal);
        if (1 > t->size[0]) {
          sigmaIsDiag = FALSE;
        } else {
          sigmaIsDiag = (t->size[0] > 2147483646);
        }

        if (sigmaIsDiag) {
          emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
          check_forloop_overflow_error();
          emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
        }

        emlrtPopRtStackR2012b(&kf_emlrtRSI, emlrtRootTLSGlobal);
        for (i = 1; i <= t->size[0]; i++) {
          bv->data[i - 1] = t->data[bk + i];
        }

        err = b_r->data[ak];
        i6 = cv->size[0];
        cv->size[0] = bv->size[0];
        emxEnsureCapacity((emxArray__common *)cv, i6, (int32_T)sizeof(real_T),
                          &fb_emlrtRTEI);
        loop_ub = bv->size[0];
        for (i6 = 0; i6 < loop_ub; i6++) {
          cv->data[i6] = err * bv->data[i6];
        }

        emlrtPushRtStackR2012b(&lf_emlrtRSI, emlrtRootTLSGlobal);
        loop_ub = t->size[0];
        if (1 > loop_ub) {
          sigmaIsDiag = FALSE;
        } else {
          loop_ub = t->size[0];
          sigmaIsDiag = (loop_ub > 2147483646);
        }

        if (sigmaIsDiag) {
          emlrtPushRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
          check_forloop_overflow_error();
          emlrtPopRtStackR2012b(&lc_emlrtRSI, emlrtRootTLSGlobal);
        }

        emlrtPopRtStackR2012b(&lf_emlrtRSI, emlrtRootTLSGlobal);
        i = 1;
        do {
          exitg2 = 0;
          loop_ub = t->size[0];
          if (i <= loop_ub) {
            r22->data[(ck + i) - 1] = cv->data[i - 1];
            i++;
          } else {
            exitg2 = 1;
          }
        } while (exitg2 == 0);

        if (asub < b_r->size[1]) {
          ak++;
          bk += t->size[0];
          bsub++;
          asub++;
        } else if (bsub < t->size[1]) {
          bk += t->size[0];
          bsub++;
        } else {
          asub = 1;
          bsub = 1;
        }

        loop_ub = t->size[0];
        ck += loop_ub;
      } else {
        exitg1 = 1;
      }
    } while (exitg1 == 0);

    emxFree_real_T(&cv);
    emxFree_real_T(&bv);
    emlrtPopRtStackR2012b(&df_emlrtRSI, emlrtRootTLSGlobal);
    for (i6 = 0; i6 < 2; i6++) {
      d_r[i6] = r22->size[i6];
    }

    for (i6 = 0; i6 < 2; i6++) {
      b_mu[i6] = mu->size[i6];
    }

    emlrtSizeEqCheck2DFastR2012b(d_r, b_mu, &cb_emlrtECI, emlrtRootTLSGlobal);
    i6 = r->size[0] * r->size[1];
    r->size[0] = r22->size[0];
    r->size[1] = r22->size[1];
    emxEnsureCapacity((emxArray__common *)r, i6, (int32_T)sizeof(real_T),
                      &fb_emlrtRTEI);
    loop_ub = r22->size[0] * r22->size[1];
    for (i6 = 0; i6 < loop_ub; i6++) {
      r->data[i6] = r22->data[i6] + mu->data[i6];
    }

    emxFree_real_T(&r22);
  } else {
    /*  Factor sigma using a function that will perform a Cholesky-like */
    /*  factorization as long as the sigma matrix is positive */
    /*  semi-definite (can have perfect correlation). Cholesky requires a */
    /*  positive definite matrix.  sigma == T'*T */
    emlrtPushRtStackR2012b(&ef_emlrtRSI, emlrtRootTLSGlobal);
    cholcov_(sigma, t, &err);
    emlrtPopRtStackR2012b(&ef_emlrtRSI, emlrtRootTLSGlobal);
    emlrtPushRtStackR2012b(&ff_emlrtRSI, emlrtRootTLSGlobal);
    emlrtPushRtStackR2012b(&gf_emlrtRSI, emlrtRootTLSGlobal);
    c_emlrt_marshallIn(randn(emlrt_marshallOut(1.0), emlrt_marshallOut(t->size[0]),
      &db_emlrtMCI), "randn", b_r);
    emlrtPopRtStackR2012b(&gf_emlrtRSI, emlrtRootTLSGlobal);
    emlrtPushRtStackR2012b(&ud_emlrtRSI, emlrtRootTLSGlobal);
    if (!(b_r->size[1] == t->size[0])) {
      if ((b_r->size[1] == 1) || ((t->size[0] == 1) && (t->size[1] == 1))) {
        emlrtPushRtStackR2012b(&wd_emlrtRSI, emlrtRootTLSGlobal);
        c_y = NULL;
        m8 = mxCreateCharArray(2, iv34);
        for (i = 0; i < 45; i++) {
          cv38[i] = cv39[i];
        }

        emlrtInitCharArrayR2013a(emlrtRootTLSGlobal, 45, m8, cv38);
        emlrtAssign(&c_y, m8);
        error(message(c_y, &r_emlrtMCI), &s_emlrtMCI);
        emlrtPopRtStackR2012b(&wd_emlrtRSI, emlrtRootTLSGlobal);
      } else {
        emlrtPushRtStackR2012b(&vd_emlrtRSI, emlrtRootTLSGlobal);
        d_y = NULL;
        m8 = mxCreateCharArray(2, iv35);
        for (i = 0; i < 21; i++) {
          cv40[i] = cv41[i];
        }

        emlrtInitCharArrayR2013a(emlrtRootTLSGlobal, 21, m8, cv40);
        emlrtAssign(&d_y, m8);
        error(message(d_y, &t_emlrtMCI), &u_emlrtMCI);
        emlrtPopRtStackR2012b(&vd_emlrtRSI, emlrtRootTLSGlobal);
      }
    }

    emlrtPopRtStackR2012b(&ud_emlrtRSI, emlrtRootTLSGlobal);
    b_emxInit_real_T(&r23, 2, &fb_emlrtRTEI, TRUE);
    if ((b_r->size[1] == 1) || (t->size[0] == 1)) {
      i6 = r->size[0] * r->size[1];
      r->size[0] = 1;
      r->size[1] = t->size[1];
      emxEnsureCapacity((emxArray__common *)r, i6, (int32_T)sizeof(real_T),
                        &fb_emlrtRTEI);
      loop_ub = t->size[1];
      for (i6 = 0; i6 < loop_ub; i6++) {
        r->data[r->size[0] * i6] = 0.0;
        i = b_r->size[1];
        for (c_r = 0; c_r < i; c_r++) {
          r->data[r->size[0] * i6] += b_r->data[b_r->size[0] * c_r] * t->
            data[c_r + t->size[0] * i6];
        }
      }
    } else {
      sz[1] = (uint32_T)t->size[1];
      i6 = r23->size[0] * r23->size[1];
      r23->size[0] = 1;
      emxEnsureCapacity((emxArray__common *)r23, i6, (int32_T)sizeof(real_T),
                        &fb_emlrtRTEI);
      i6 = r23->size[0] * r23->size[1];
      r23->size[1] = (int32_T)sz[1];
      emxEnsureCapacity((emxArray__common *)r23, i6, (int32_T)sizeof(real_T),
                        &fb_emlrtRTEI);
      loop_ub = (int32_T)sz[1];
      for (i6 = 0; i6 < loop_ub; i6++) {
        r23->data[i6] = 0.0;
      }

      emlrtPushRtStackR2012b(&ge_emlrtRSI, emlrtRootTLSGlobal);
      emlrtPushRtStackR2012b(&he_emlrtRSI, emlrtRootTLSGlobal);
      i6 = r->size[0] * r->size[1];
      r->size[0] = 1;
      emxEnsureCapacity((emxArray__common *)r, i6, (int32_T)sizeof(real_T),
                        &fb_emlrtRTEI);
      i6 = r->size[0] * r->size[1];
      r->size[1] = (int32_T)sz[1];
      emxEnsureCapacity((emxArray__common *)r, i6, (int32_T)sizeof(real_T),
                        &fb_emlrtRTEI);
      loop_ub = (int32_T)sz[1];
      for (i6 = 0; i6 < loop_ub; i6++) {
        r->data[i6] = 0.0;
      }

      if ((t->size[1] < 1) || (b_r->size[1] < 1)) {
      } else {
        emlrtPushRtStackR2012b(&je_emlrtRSI, emlrtRootTLSGlobal);
        err = 1.0;
        beta1 = 0.0;
        TRANSB = 'N';
        TRANSA = 'N';
        emlrtPushRtStackR2012b(&ne_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        m_t = (ptrdiff_t)(1);
        emlrtPopRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPopRtStackR2012b(&ne_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&oe_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        n_t = (ptrdiff_t)(t->size[1]);
        emlrtPopRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPopRtStackR2012b(&oe_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&pe_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        k_t = (ptrdiff_t)(b_r->size[1]);
        emlrtPopRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPopRtStackR2012b(&pe_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&qe_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        lda_t = (ptrdiff_t)(1);
        emlrtPopRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPopRtStackR2012b(&qe_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&re_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        ldb_t = (ptrdiff_t)(b_r->size[1]);
        emlrtPopRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPopRtStackR2012b(&re_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&se_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        ldc_t = (ptrdiff_t)(1);
        emlrtPopRtStackR2012b(&af_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPopRtStackR2012b(&se_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&te_emlrtRSI, emlrtRootTLSGlobal);
        alpha1_t = (double *)(&err);
        emlrtPopRtStackR2012b(&te_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&ue_emlrtRSI, emlrtRootTLSGlobal);
        Aia0_t = (double *)(&b_r->data[0]);
        emlrtPopRtStackR2012b(&ue_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&ve_emlrtRSI, emlrtRootTLSGlobal);
        Bib0_t = (double *)(&t->data[0]);
        emlrtPopRtStackR2012b(&ve_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&we_emlrtRSI, emlrtRootTLSGlobal);
        beta1_t = (double *)(&beta1);
        emlrtPopRtStackR2012b(&we_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&xe_emlrtRSI, emlrtRootTLSGlobal);
        Cic0_t = (double *)(&r->data[0]);
        emlrtPopRtStackR2012b(&xe_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPushRtStackR2012b(&ye_emlrtRSI, emlrtRootTLSGlobal);
        dgemm(&TRANSA, &TRANSB, &m_t, &n_t, &k_t, alpha1_t, Aia0_t, &lda_t,
              Bib0_t, &ldb_t, beta1_t, Cic0_t, &ldc_t);
        emlrtPopRtStackR2012b(&ye_emlrtRSI, emlrtRootTLSGlobal);
        emlrtPopRtStackR2012b(&je_emlrtRSI, emlrtRootTLSGlobal);
      }

      emlrtPopRtStackR2012b(&he_emlrtRSI, emlrtRootTLSGlobal);
      emlrtPopRtStackR2012b(&ge_emlrtRSI, emlrtRootTLSGlobal);
    }

    emxFree_real_T(&r23);
    emlrtPopRtStackR2012b(&ff_emlrtRSI, emlrtRootTLSGlobal);
    for (i6 = 0; i6 < 2; i6++) {
      d_r[i6] = r->size[i6];
    }

    for (i6 = 0; i6 < 2; i6++) {
      b_mu[i6] = mu->size[i6];
    }

    emlrtSizeEqCheck2DFastR2012b(d_r, b_mu, &db_emlrtECI, emlrtRootTLSGlobal);
    i6 = r->size[0] * r->size[1];
    r->size[0] = 1;
    emxEnsureCapacity((emxArray__common *)r, i6, (int32_T)sizeof(real_T),
                      &fb_emlrtRTEI);
    i = r->size[0];
    c_r = r->size[1];
    loop_ub = i * c_r;
    for (i6 = 0; i6 < loop_ub; i6++) {
      r->data[i6] += mu->data[i6];
    }
  }

  emxFree_real_T(&b_r);
  emxFree_real_T(&t);

  /*  Multiple covariance matrices */
  emlrtHeapReferenceStackLeaveFcnR2012b(emlrtRootTLSGlobal);
}

/* End of code generation (mvnrnd_.c) */