aa.c in scs-0.3.2

- old
+ new

@@ -94,10 +94,14 @@
 #define TIME_TIC
 #define TIME_TOC
 
 #endif
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 /* BLAS functions used */
 aa_float BLAS(nrm2)(blas_int *n, aa_float *x, blas_int *incx);
 void BLAS(axpy)(blas_int *n, aa_float *a, const aa_float *x, blas_int *incx,
                 aa_float *y, blas_int *incy);
 void BLAS(gemv)(const char *trans, const blas_int *m, const blas_int *n,
@@ -111,10 +115,14 @@
                 blas_int *lda, aa_float *b, blas_int *ldb, aa_float *beta,
                 aa_float *c, blas_int *ldc);
 void BLAS(scal)(const blas_int *n, const aa_float *a, aa_float *x,
                 const blas_int *incx);
 
+#ifdef __cplusplus
+}
+#endif
+
 /* This file uses Anderson acceleration to improve the convergence of
  * a fixed point mapping.
  * At each iteration we need to solve a (small) linear system, we
  * do this using LAPACK ?gesv.
  */
@@ -274,21 +282,22 @@
 }
 
 /* f = (1-relaxation) * \sum_i a_i x_i + relaxation * \sum_i a_i f_i */
 static void relax(aa_float *f, AaWork *a, aa_int len) {
   TIME_TIC
-  /* x_work = x - S * work */
+  /* x_work = x initially */
   blas_int bdim = (blas_int)(a->dim), one = 1, blen = (blas_int)len;
   aa_float onef = 1.0, neg_onef = -1.0;
   aa_float one_m_relaxation = 1. - a->relaxation;
+  /* x_work = x - S * work */
   BLAS(gemv)
   ("NoTrans", &bdim, &blen, &neg_onef, a->S, &bdim, a->work, &one, &onef,
    a->x_work, &one);
   /* f = relaxation * f */
-  BLAS(scal)(&blen, &a->relaxation, f, &one);
+  BLAS(scal)(&bdim, &a->relaxation, f, &one);
   /* f += (1 - relaxation) * x_work */
-  BLAS(axpy)(&blen, &one_m_relaxation, a->x_work, &one, f, &one);
+  BLAS(axpy)(&bdim, &one_m_relaxation, a->x_work, &one, f, &one);
   TIME_TOC
 }
 
 /* solves the system of equations to perform the AA update
  * at the end f contains the next iterate to be returned
@@ -350,10 +359,10 @@
                 aa_float max_weight_norm, aa_int verbosity) {
   TIME_TIC
   AaWork *a = (AaWork *)calloc(1, sizeof(AaWork));
   if (!a) {
     printf("Failed to allocate memory for AA.\n");
-    return (void *)0;
+    return (AaWork *)0;
   }
   a->type1 = type1;
   a->iter = 0;
   a->dim = dim;
   a->mem = MIN(mem, dim); /* for rank stability */