/* [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "meq.h"

# ifdef SINGLE
# define FCT sbslv
# define BCP cblas_scopy
# define BMV cblas_sgemv
# else
# define BCP cblas_dcopy
# define BMV cblas_dgemv
# define FCT dbslv
# endif
# define CPY(y,x,n) BCP(n,x,1,y,1) /* y = x */
# define MVM(y,A,x,m,n) BMV(CblasColMajor, CblasNoTrans, m, n, FLTC(-1.0), A, m, x, 1, FLTC(1.0), y, 1) /* y = -A*x + y */

void FCT(FLT *a, FLT *b, FLT *a1, FLT *b1, FLT *b0, FLT *A, FLT *B, int na, int nb, int niter) {
 int kiter = niter;
 /* b = b0 */
 CPY(b,b0,nb);
 for (;kiter--;) {
  /* a = a1 - A*b */
  CPY(a,a1,na); /* a = a1 */
  MVM(a,A,b,na,nb); /* a = a - A*b */
  /* b = b1 -B*a */
  CPY(b,b1,nb); /* b = b1 */
  MVM(b,B,a,nb,na); /* b = b - B*a */
 }
}
