text: #include <ulmblas.h>
  #include <auxiliary/xerbla.h>
  #include <math.h>
  
  void
  ULMBLAS(dsymm)(const enum Side   side,
                 const enum UpLo   upLo,
                 const int         m,
                 const int         n,
                 const double      alpha,
                 const double      *A,
                 const int         ldA,
                 const double      *B,
                 const int         ldB,
                 const double      beta,
                 double            *C,
                 const int         ldC)
  {
  //
  //  Local scalars
  //
      int     i, j, k;
      double  tmp1, tmp2;
  
  //
  //  Quick return if possible.
  //
      if (m==0 || n==0 || (alpha==0.0 && beta==1.0)) {
          return;
      }
  
  //
  //  And when alpha is exactly zero
  //
      if (alpha==0.0) {
          if (beta==0.0) {
              for (j=0; j<n; ++j) {
                  for (i=0; i<m; ++i) {
                      C[i+j*ldC] = 0.0;
                  }
              }
          } else {
              for (j=0; j<n; ++j) {
                  for (i=0; i<m; ++i) {
                      C[i+j*ldC] *= beta;
                  }
              }
          }
          return;
      }
  //
  //  Start the operations.
  //
      if (side==Left) {
  //
  //      Form  C := alpha*A*B + beta*C.
  //
          if (upLo==Upper) {
              for (j=0; j<n; ++j) {
                  for (i=0; i<m; ++i) {
                      tmp1 = alpha*B[i+j*ldB];
                      tmp2 = 0.0;
                      for (k=0; k<i; ++k) {
                          C[k+j*ldC] += tmp1*A[k+i*ldA];
                          tmp2       += B[k+j*ldB]*A[k+i*ldA];
                      }
                      if (beta==0.0) {
                          C[i+j*ldC] = tmp1*A[i+i*ldA] + alpha*tmp2;
                      } else {
                          C[i+j*ldC] = beta*C[i+j*ldC]
                                     + tmp1*A[i+i*ldA]
                                     + alpha*tmp2;
                      }
                  }
              }
          } else {
              for (j=0; j<n; ++j) {
                  for (i=m-1; i>=0; --i) {
                      tmp1 = alpha*B[i+j*ldB];
                      tmp2 = 0.0;
                      for (k=i+1; k<m; ++k) {
                          C[k+j*ldC] += tmp1*A[k+i*ldA];
                          tmp2       += B[k+j*ldB]*A[k+i*ldA];
                      }
                      if (beta==0.0) {
                          C[i+j*ldC] = tmp1*A[i+i*ldA] + alpha*tmp2;
                      } else {
                          C[i+j*ldC] = beta*C[i+j*ldC]
                                     + tmp1*A[i+i*ldA]
                                     + alpha*tmp2;
                      }
                   }
              }
          }
      } else {
  //
  //      Form  C := alpha*B*A + beta*C.
  //
          for (j=0; j<n; ++j) {
              tmp1 = alpha*A[j+j*ldA];
              if (beta==0.0) {
                  for (i=0; i<m; ++i) {
                      C[i+j*ldC] = tmp1*B[i+j*ldB];
                  }
              } else {
                  for (i=0; i<m; ++i) {
                      C[i+j*ldC] = beta*C[i+j*ldC] + tmp1*B[i+j*ldB];
                  }
              }
              for (k=0; k<j; ++k) {
                  if (upLo==Upper) {
                      tmp1 = alpha*A[k+j*ldA];
                  } else {
                      tmp1 = alpha*A[j+k*ldA];
                  }
                  for (i=0; i<m; ++i) {
                      C[i+j*ldC] += tmp1*B[i+k*ldB];
                  }
              }
              for (k=j+1; k<n; ++k) {
                  if (upLo==Upper) {
                      tmp1 = alpha*A[j+k*ldA];
                  } else {
                      tmp1 = alpha*A[k+j*ldA];
                  }
                  for (i=0; i<m; ++i) {
                      C[i+j*ldC] += tmp1*B[i+k*ldB];
                  }
              }
          }
      }
  }
  
  void
  F77BLAS(dsymm)(const char     *_side,
                 const char     *_upLo,
                 const int      *_m,
                 const int      *_n,
                 const double   *_alpha,
                 const double   *A,
                 const int      *_ldA,
                 const double   *B,
                 const int      *_ldB,
                 const double   *_beta,
                 double         *C,
                 const int      *_ldC)
  {
  //
  //  Dereference scalar parameters
  //
      enum Side side = charToSide(*_side);
      enum UpLo upLo = charToUpLo(*_upLo);
      int m          = *_m;
      int n          = *_n;
      double alpha   = *_alpha;
      int ldA        = *_ldA;
      int ldB        = *_ldB;
      double beta    = *_beta;
      int ldC        = *_ldC;
  
  //
  //  Set  numRowsA as the number of rows of A
  //
      int numRowsA = (side==Left) ? m : n;
  
  //
  //  Test the input parameters
  //
      int info = 0;
      if (side==0) {
          info = 1;
      } else if (upLo==0) {
          info = 2;
      } else if (m<0) {
          info = 3;
      } else if (n<0) {
          info = 4;
      } else if (ldA<max(1,numRowsA)) {
          info = 7;
      } else if (ldB<m) {
          info = 9;
      } else if (ldC<m) {
          info = 12;
      }
  
      if (info!=0) {
          F77BLAS(xerbla)("DSYMM ", &info);
      }
  
      ULMBLAS(dsymm)(side, upLo,
                     m, n,
                     alpha,
                     A, ldA,
                     B, ldB,
                     beta,
                     C, ldC);
  }
  