cblas_base.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2018, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /*
24  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
25  * Copyright (C) 2001, 2007 Brian Gough
26  *
27  * This program is free software; you can redistribute it and/or modify
28  * it under the terms of the GNU General Public License as published by
29  * the Free Software Foundation; either version 3 of the License, or (at
30  * your option) any later version.
31  *
32  * This program is distributed in the hope that it will be useful, but
33  * WITHOUT ANY WARRANTY; without even the implied warranty of
34  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
35  * General Public License for more details.
36  *
37  * You should have received a copy of the GNU General Public License
38  * along with this program; if not, write to the Free Software
39  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
40  * 02110-1301, USA.
41  */
42 /** \file cblas_base.h
43  \brief O2scl basic linear algebra function templates
44 
45  See \ref o2scl_cblas for more documentation on
46  these functions.
47 
48  \future Add float and complex versions?
49  \future There are some Level-1 BLAS functions which are already
50  present in vector.h. Should we move all of them to one place or
51  the other? The ones in vector.h are generic in the sense that they
52  can use doubles or floats, but the ones here can use either () or
53  [].
54 
55 */
56 
57 #ifdef DOXYGEN
58 namespace o2scl_cblas {
59 #endif
60 
61  /// Matrix order, either column-major or row-major
62  enum o2cblas_order {o2cblas_RowMajor=101, o2cblas_ColMajor=102};
63 
64  /// Transpose operations
65  enum o2cblas_transpose {o2cblas_NoTrans=111, o2cblas_Trans=112,
66  o2cblas_ConjTrans=113};
67 
68  /// Upper- or lower-triangular
69  enum o2cblas_uplo {o2cblas_Upper=121, o2cblas_Lower=122};
70 
71  /// Unit or generic diagonal
72  enum o2cblas_diag {o2cblas_NonUnit=131, o2cblas_Unit=132};
73 
74  /// Left or right sided operation
75  enum o2cblas_side {o2cblas_Left=141, o2cblas_Right=142};
76 
77  /// \name Test if matrix elements are finite
78  //@{
79  /** \brief Test if the first \c n elements of a matrix are finite
80 
81  If \c n is zero, this will return true without throwing
82  an exception.
83  */
84  template<class mat_t>
85  bool matrix_is_finite(size_t m, size_t n, mat_t &data) {
86  for(size_t i=0;i<m;i++) {
87  for(size_t j=0;j<n;j++) {
88  if (!std::isfinite(O2SCL_IX2(data,i,j))) return false;
89  }
90  }
91  return true;
92  }
93 
94  /** \brief Test if a matrix is finite
95 
96  If \c n is zero, this will return true without throwing
97  an exception.
98  */
99  template<class mat_t> bool matrix_is_finite(mat_t &data) {
100  return matrix_is_finite(data.size1(),data.size2(),data);
101  }
102  //@}
103 
104  /// \name Level-1 BLAS functions
105  //@{
106  /** \brief Compute the absolute sum of vector elements
107 
108  If \c alpha is zero, this function returns and performs
109  no computations.
110  */
111  template<class vec_t> double dasum(const size_t N, const vec_t &X) {
112  double r=0.0;
113  for(size_t i=0;i<N;i++) {
114  r+=fabs(X[i]);
115  }
116  return r;
117  }
118 
119  /** \brief Compute \f$ y= \alpha x+y \f$
120 
121  If \c alpha is zero, this function returns and performs
122  no computations.
123  */
124  template<class vec_t, class vec2_t>
125  void daxpy(const double alpha, const size_t N, const vec_t &X,
126  vec2_t &Y) {
127 
128  size_t i;
129 
130  if (alpha == 0.0) {
131  return;
132  }
133 
134  const size_t m=N % 4;
135 
136  for (i=0;i<m;i++) {
137  O2SCL_IX(Y,i)+=alpha*O2SCL_IX(X,i);
138  }
139 
140  for (i=m;i+3<N;i+=4) {
141  O2SCL_IX(Y,i)+=alpha*O2SCL_IX(X,i);
142  O2SCL_IX(Y,i+1)+=alpha*O2SCL_IX(X,i+1);
143  O2SCL_IX(Y,i+2)+=alpha*O2SCL_IX(X,i+2);
144  O2SCL_IX(Y,i+3)+=alpha*O2SCL_IX(X,i+3);
145  }
146  }
147 
148  /// Compute \f$ r=x \cdot y \f$
149  template<class vec_t, class vec2_t>
150  double ddot(const size_t N, const vec_t &X, const vec2_t &Y) {
151 
152  double r=0.0;
153  size_t i;
154 
155  const size_t m=N % 4;
156 
157  for (i=0;i<m;i++) {
158  r+=O2SCL_IX(X,i)*O2SCL_IX(Y,i);
159  }
160 
161  for (i=m;i+3<N;i+=4) {
162  r+=O2SCL_IX(X,i)*O2SCL_IX(Y,i);
163  r+=O2SCL_IX(X,i+1)*O2SCL_IX(Y,i+1);
164  r+=O2SCL_IX(X,i+2)*O2SCL_IX(Y,i+2);
165  r+=O2SCL_IX(X,i+3)*O2SCL_IX(Y,i+3);
166  }
167 
168  return r;
169  }
170 
171  /** \brief Compute the norm of the vector \c X
172 
173  \note The suffix "2" on the function name indicates that this
174  computes the "2-norm", not that the norm is squared.
175 
176  If \c N is less than or equal to zero, this function returns
177  zero without calling the error handler.
178 
179  This function works only with vectors which hold \c double. For
180  the norm of a general floating point vector, see \ref
181  vector_norm().
182  */
183  template<class vec_t> double dnrm2(const size_t N, const vec_t &X) {
184 
185  double scale=0.0;
186  double ssq=1.0;
187  size_t i;
188 
189  if (N == 0) {
190  return 0;
191  } else if (N == 1) {
192  return fabs(O2SCL_IX(X,0));
193  }
194 
195  for (i=0;i<N;i++) {
196  const double x=O2SCL_IX(X,i);
197 
198  if (x != 0.0) {
199  const double ax=fabs(x);
200 
201  if (scale<ax) {
202  ssq=1.0+ssq*(scale/ax)*(scale/ax);
203  scale=ax;
204  } else {
205  ssq+=(ax/scale)*(ax/scale);
206  }
207  }
208 
209  }
210 
211  return scale*sqrt(ssq);
212  }
213 
214  /** \brief Compute \f$ x=\alpha x \f$
215  */
216  template<class vec_t>
217  void dscal(const double alpha, const size_t N, vec_t &X) {
218 
219  size_t i;
220  const size_t m=N % 4;
221 
222  for (i=0;i<m;i++) {
223  O2SCL_IX(X,i)*=alpha;
224  }
225 
226  for (i=m;i+3<N;i+=4) {
227  O2SCL_IX(X,i)*=alpha;
228  O2SCL_IX(X,i+1)*=alpha;
229  O2SCL_IX(X,i+2)*=alpha;
230  O2SCL_IX(X,i+3)*=alpha;
231  }
232  }
233  //@}
234 
235  /// \name Level-2 BLAS functions
236  //@{
237  /** \brief Compute \f$ y=\alpha \left[\mathrm{op}(A)\right] x+
238  \beta y \f$.
239 
240  If \c M or \c N is zero, or if \c alpha is zero and \c beta is
241  one, this function performs no calculations and returns without
242  calling the error handler.
243  */
244  template<class mat_t, class vec_t, class vec2_t>
245  void dgemv(const enum o2cblas_order order,
246  const enum o2cblas_transpose TransA, const size_t M,
247  const size_t N, const double alpha, const mat_t &A,
248  const vec_t &X, const double beta, vec2_t &Y) {
249 
250  size_t i, j;
251  size_t lenX, lenY;
252 
253  // If conjugate transpose is requested, just assume plain transpose
254  const int Trans=(TransA != o2cblas_ConjTrans) ? TransA : o2cblas_Trans;
255 
256  if (M == 0 || N == 0) {
257  return;
258  }
259 
260  if (alpha == 0.0 && beta == 1.0) {
261  return;
262  }
263 
264  if (Trans == o2cblas_NoTrans) {
265  lenX=N;
266  lenY=M;
267  } else {
268  lenX=M;
269  lenY=N;
270  }
271 
272  /* form y := beta*y */
273  if (beta == 0.0) {
274  size_t iy=0;
275  for (i=0;i<lenY;i++) {
276  O2SCL_IX(Y,iy)=0.0;
277  iy++;
278  }
279  } else if (beta != 1.0) {
280  size_t iy=0;
281  for (i=0;i<lenY;i++) {
282  O2SCL_IX(Y,iy) *= beta;
283  iy++;
284  }
285  }
286 
287  if (alpha == 0.0) {
288  return;
289  }
290 
291  if ((order == o2cblas_RowMajor && Trans == o2cblas_NoTrans) ||
292  (order == o2cblas_ColMajor && Trans == o2cblas_Trans)) {
293 
294  /* form y := alpha*A*x+y */
295  size_t iy=0;
296  for (i=0;i<lenY;i++) {
297  double temp=0.0;
298  size_t ix=0;
299  for (j=0;j<lenX;j++) {
300  temp+=O2SCL_IX(X,ix)*O2SCL_IX2(A,i,j);
301  ix++;
302  }
303  O2SCL_IX(Y,iy)+=alpha*temp;
304  iy++;
305  }
306 
307  } else if ((order == o2cblas_RowMajor && Trans == o2cblas_Trans) ||
308  (order == o2cblas_ColMajor && Trans == o2cblas_NoTrans)) {
309 
310  /* form y := alpha*A'*x+y */
311  size_t ix=0;
312  for (j=0;j<lenX;j++) {
313  const double temp=alpha*O2SCL_IX(X,ix);
314  if (temp != 0.0) {
315  size_t iy=0;
316  for (i=0;i<lenY;i++) {
317  O2SCL_IX(Y,iy)+=temp*O2SCL_IX2(A,j,i);
318  iy++;
319  }
320  }
321  ix++;
322  }
323 
324  } else {
325  O2SCL_ERR("Unrecognized operation in dgemv().",o2scl::exc_einval);
326  }
327  return;
328  }
329 
330  /** \brief Compute \f$ x=\mathrm{op} (A)^{-1} x \f$
331 
332  If \c N is zero, this function does nothing and returns zero.
333  */
334  template<class mat_t, class vec_t>
335  void dtrsv(const enum o2cblas_order order,
336  const enum o2cblas_uplo Uplo,
337  const enum o2cblas_transpose TransA,
338  const enum o2cblas_diag Diag,
339  const size_t M, const size_t N, const mat_t &A, vec_t &X) {
340 
341  const int nonunit=(Diag == o2cblas_NonUnit);
342  int ix, jx;
343  int i, j;
344  const int Trans=(TransA != o2cblas_ConjTrans) ? TransA : o2cblas_Trans;
345 
346  if (N == 0) {
347  return;
348  }
349 
350  /* form x := inv( A )*x */
351 
352  if ((order == o2cblas_RowMajor && Trans == o2cblas_NoTrans &&
353  Uplo == o2cblas_Upper) ||
354  (order == o2cblas_ColMajor && Trans == o2cblas_Trans &&
355  Uplo == o2cblas_Lower)) {
356 
357  /* backsubstitution */
358 
359  // O2scl: Note that subtraction of 1 from size_t N here is ok
360  // since we already handled the N=0 case above
361  ix=(int)(N-1);
362  if (nonunit) {
363  O2SCL_IX(X,ix)=O2SCL_IX(X,ix)/O2SCL_IX2(A,N-1,N-1);
364  }
365  ix--;
366 
367  for (i=(int)(N-1);i>0 && i--;) {
368  double tmp=O2SCL_IX(X,ix);
369  jx=ix +1;
370  for (j=i+1;j<((int)N);j++) {
371  const double Aij=O2SCL_IX2(A,i,j);
372  tmp-=Aij*O2SCL_IX(X,jx);
373  jx++;
374  }
375  if (nonunit) {
376  O2SCL_IX(X,ix)=tmp/O2SCL_IX2(A,i,i);
377  } else {
378  O2SCL_IX(X,ix)=tmp;
379  }
380  ix--;
381  }
382 
383  } else if ((order == o2cblas_RowMajor && Trans == o2cblas_NoTrans &&
384  Uplo == o2cblas_Lower) ||
385  (order == o2cblas_ColMajor && Trans == o2cblas_Trans &&
386  Uplo == o2cblas_Upper)) {
387 
388  /* forward substitution */
389  ix=0;
390  if (nonunit) {
391  O2SCL_IX(X,ix)=O2SCL_IX(X,ix)/O2SCL_IX2(A,0,0);
392  }
393  ix++;
394  for (i=1;i<((int)N);i++) {
395  double tmp=O2SCL_IX(X,ix);
396  jx=0;
397  for (j=0;j<i;j++) {
398  const double Aij=O2SCL_IX2(A,i,j);
399  tmp-=Aij*O2SCL_IX(X,jx);
400  jx++;
401  }
402  if (nonunit) {
403  O2SCL_IX(X,ix)=tmp/O2SCL_IX2(A,i,i);
404  } else {
405  O2SCL_IX(X,ix)=tmp;
406  }
407  ix++;
408  }
409 
410  } else if ((order == o2cblas_RowMajor && Trans == o2cblas_Trans &&
411  Uplo == o2cblas_Upper) ||
412  (order == o2cblas_ColMajor && Trans == o2cblas_NoTrans &&
413  Uplo == o2cblas_Lower)) {
414 
415  /* form x := inv( A' )*x */
416 
417  /* forward substitution */
418  ix=0;
419  if (nonunit) {
420  O2SCL_IX(X,ix)=O2SCL_IX(X,ix)/O2SCL_IX2(A,0,0);
421  }
422  ix++;
423  for (i=1;i<((int)N);i++) {
424  double tmp=O2SCL_IX(X,ix);
425  jx=0;
426  for (j=0;j<i;j++) {
427  const double Aji=O2SCL_IX2(A,j,i);
428  tmp-=Aji*O2SCL_IX(X,jx);
429  jx++;
430  }
431  if (nonunit) {
432  O2SCL_IX(X,ix)=tmp/O2SCL_IX2(A,i,i);
433  } else {
434  O2SCL_IX(X,ix)=tmp;
435  }
436  ix++;
437  }
438 
439  } else if ((order == o2cblas_RowMajor && Trans == o2cblas_Trans &&
440  Uplo == o2cblas_Lower) ||
441  (order == o2cblas_ColMajor && Trans == o2cblas_NoTrans &&
442  Uplo == o2cblas_Upper)) {
443 
444  /* backsubstitution */
445  // O2scl: Note that subtraction of 1 from size_t N here is ok
446  // since we already handled the N=0 case above
447  ix=(int)(N-1);
448  if (nonunit) {
449  O2SCL_IX(X,ix)=O2SCL_IX(X,ix)/O2SCL_IX2(A,N-1,N-1);
450  }
451  ix--;
452  for (i=(int)(N-1);i>0 && i--;) {
453  double tmp=O2SCL_IX(X,ix);
454  jx=ix+1;
455  for (j=i+1;j<((int)N);j++) {
456  const double Aji=O2SCL_IX2(A,j,i);
457  tmp-=Aji*O2SCL_IX(X,jx);
458  jx++;
459  }
460  if (nonunit) {
461  O2SCL_IX(X,ix)=tmp/O2SCL_IX2(A,i,i);
462  } else {
463  O2SCL_IX(X,ix)=tmp;
464  }
465  ix--;
466  }
467 
468  } else {
469  O2SCL_ERR("Unrecognized operation in dtrsv().",o2scl::exc_einval);
470  }
471  return;
472  }
473 
474  /** \brief Compute \f$ x=op(A) x \f$ for the triangular matrix \c A
475  */
476  template<class mat_t, class vec_t>
477  void dtrmv(const enum o2cblas_order Order,
478  const enum o2cblas_uplo Uplo,
479  const enum o2cblas_transpose TransA,
480  const enum o2cblas_diag Diag, const size_t N,
481  const mat_t &A, vec_t &x) {
482 
483  int i, j;
484 
485  const int nonunit=(Diag == o2cblas_NonUnit);
486  const int Trans=(TransA != o2cblas_ConjTrans) ? TransA : o2cblas_Trans;
487 
488  if ((Order == o2cblas_RowMajor && Trans == o2cblas_NoTrans &&
489  Uplo == o2cblas_Upper) ||
490  (Order == o2cblas_ColMajor && Trans == o2cblas_Trans &&
491  Uplo == o2cblas_Lower)) {
492 
493  /* form x := A*x */
494 
495  for (i=0;i<N;i++) {
496  double temp=0.0;
497  const size_t j_min=i+1;
498  const size_t j_max=N;
499  size_t jx=j_min;
500  for (j=j_min;j<j_max;j++) {
501  temp+=O2SCL_IX(x,jx)*O2SCL_IX2(A,i,j);
502  jx++;
503  }
504  if (nonunit) {
505  O2SCL_IX(x,i)=temp+O2SCL_IX(x,i)*O2SCL_IX2(A,i,i);
506  } else {
507  O2SCL_IX(x,i)+=temp;
508  }
509  }
510 
511  } else if ((Order == o2cblas_RowMajor && Trans == o2cblas_NoTrans &&
512  Uplo == o2cblas_Lower) ||
513  (Order == o2cblas_ColMajor && Trans == o2cblas_Trans &&
514  Uplo == o2cblas_Upper)) {
515 
516  // O2scl: Note that subtraction of 1 from size_t N here is ok
517  // since we already handled the N=0 case above
518  size_t ix=N-1;
519  for (i=N;i>0 && i--;) {
520  double temp=0.0;
521  const size_t j_min=0;
522  const size_t j_max=i;
523  size_t jx=j_min;
524  for (j=j_min;j<j_max;j++) {
525  temp+=O2SCL_IX(x,jx)*O2SCL_IX2(A,i,j);
526  jx++;
527  }
528  if (nonunit) {
529  O2SCL_IX(x,ix)=temp+O2SCL_IX(x,ix)*O2SCL_IX2(A,i,i);
530  } else {
531  O2SCL_IX(x,ix)+=temp;
532  }
533  ix--;
534  }
535 
536  } else if ((Order == o2cblas_RowMajor && Trans == o2cblas_Trans &&
537  Uplo == o2cblas_Upper) ||
538  (Order == o2cblas_ColMajor && Trans == o2cblas_NoTrans &&
539  Uplo == o2cblas_Lower)) {
540 
541  /* form x := A'*x */
542  size_t ix=N-1;
543  for (i=N;i>0 && i--;) {
544  double temp=0.0;
545  const size_t j_min=0;
546  const size_t j_max=i;
547  size_t jx=j_min;
548  for (j=j_min;j<j_max;j++) {
549  temp+=O2SCL_IX(x,jx)*O2SCL_IX2(A,j,i);
550  jx++;
551  }
552  if (nonunit) {
553  O2SCL_IX(x,ix)=temp+O2SCL_IX(x,ix)*O2SCL_IX2(A,i,i);
554  } else {
555  O2SCL_IX(x,ix)+=temp;
556  }
557  ix--;
558  }
559 
560  } else if ((Order == o2cblas_RowMajor && Trans == o2cblas_Trans &&
561  Uplo == o2cblas_Lower) ||
562  (Order == o2cblas_ColMajor && Trans == o2cblas_NoTrans &&
563  Uplo == o2cblas_Upper)) {
564 
565  for (i=0;i<N;i++) {
566  double temp=0.0;
567  const size_t j_min=i+1;
568  const size_t j_max=N;
569  size_t jx=i+1;
570  for (j=j_min;j<j_max;j++) {
571  temp+=O2SCL_IX(x,jx)*O2SCL_IX2(A,j,i);
572  jx++;
573  }
574  if (nonunit) {
575  O2SCL_IX(x,i)=temp+O2SCL_IX(x,i)*O2SCL_IX2(A,i,i);
576  } else {
577  O2SCL_IX(x,i)+=temp;
578  }
579  }
580 
581  } else {
582  O2SCL_ERR("Unrecognized operation in dtrmv().",
584  }
585 
586  return;
587  }
588  //@}
589 
590 
591  /// \name Level-3 BLAS functions
592  //@{
593  /** \brief Compute \f$ y=\alpha \mathrm{op}(A) \mathrm{op}(B) +
594  \beta C \f$
595 
596  \comment
597  If \c Order is \c RowMajor, then the matrix \c A
598  has \c M rows and \c K columns, the matrix \c B has
599  \c K rows and \c N columns, and the
600  matrix \c C has \c M rows and \c N columns.
601 
602  Is this right?
603  \endcommment
604 
605  This function works for all values of \c Order, \c TransA, and
606  \c TransB.
607  */
608  template<class mat_t>
609  void dgemm(const enum o2cblas_order Order,
610  const enum o2cblas_transpose TransA,
611  const enum o2cblas_transpose TransB, const size_t M,
612  const size_t N, const size_t K, const double alpha,
613  const mat_t &A, const mat_t &B, const double beta, mat_t &C) {
614 
615  size_t i, j, k;
616  size_t n1, n2;
617  int TransF, TransG;
618 
619  if (alpha == 0.0 && beta == 1.0) {
620  return;
621  }
622 
623  /*
624  This is a little more complicated than the original in GSL,
625  which assigned the matrices A and B to variables *F and *G which
626  then allowed some code duplication. We can't really do that
627  here, since we don't have that kind of type info on A and B, so
628  we just handle the two cases separately.
629  */
630 
631  if (Order == o2cblas_RowMajor) {
632 
633  n1=M;
634  n2=N;
635 
636  /* form y := beta*y */
637  if (beta == 0.0) {
638  for (i=0;i<n1;i++) {
639  for (j=0;j<n2;j++) {
640  O2SCL_IX2(C,i,j)=0.0;
641  }
642  }
643  } else if (beta != 1.0) {
644  for (i=0;i<n1;i++) {
645  for (j=0;j<n2;j++) {
646  O2SCL_IX2(C,i,j)*=beta;
647  }
648  }
649  }
650 
651  if (alpha == 0.0) {
652  return;
653  }
654 
655  TransF=(TransA == o2cblas_ConjTrans) ? o2cblas_Trans : TransA;
656  TransG=(TransB == o2cblas_ConjTrans) ? o2cblas_Trans : TransB;
657 
658  if (TransF == o2cblas_NoTrans && TransG == o2cblas_NoTrans) {
659 
660  /* form C := alpha*A*B+C */
661 
662  for (k=0;k<K;k++) {
663  for (i=0;i<n1;i++) {
664  const double temp=alpha*O2SCL_IX2(A,i,k);
665  if (temp != 0.0) {
666  for (j=0;j<n2;j++) {
667  O2SCL_IX2(C,i,j)+=temp*O2SCL_IX2(B,k,j);
668  }
669  }
670  }
671  }
672 
673  } else if (TransF == o2cblas_NoTrans && TransG == o2cblas_Trans) {
674 
675  /* form C := alpha*A*B'+C */
676 
677  for (i=0;i<n1;i++) {
678  for (j=0;j<n2;j++) {
679  double temp=0.0;
680  for (k=0;k<K;k++) {
681  temp+=O2SCL_IX2(A,i,k)*O2SCL_IX2(B,j,k);
682  }
683  O2SCL_IX2(C,i,j)+=alpha*temp;
684  }
685  }
686 
687  } else if (TransF == o2cblas_Trans && TransG == o2cblas_NoTrans) {
688 
689  for (k=0;k<K;k++) {
690  for (i=0;i<n1;i++) {
691  const double temp=alpha*O2SCL_IX2(A,k,i);
692  if (temp != 0.0) {
693  for (j=0;j<n2;j++) {
694  O2SCL_IX2(C,i,j)+=temp*O2SCL_IX2(B,k,j);
695  }
696  }
697  }
698  }
699 
700  } else if (TransF == o2cblas_Trans && TransG == o2cblas_Trans) {
701 
702  for (i=0;i<n1;i++) {
703  for (j=0;j<n2;j++) {
704  double temp=0.0;
705  for (k=0;k<K;k++) {
706  temp+=O2SCL_IX2(A,k,i)*O2SCL_IX2(B,j,k);
707  }
708  O2SCL_IX2(C,i,j)+=alpha*temp;
709  }
710  }
711 
712  } else {
713  O2SCL_ERR("Unrecognized operation in dgemm().",o2scl::exc_einval);
714  }
715 
716  } else {
717 
718  // Column-major case
719 
720  n1=N;
721  n2=M;
722 
723  /* form y := beta*y */
724  if (beta == 0.0) {
725  for (i=0;i<n1;i++) {
726  for (j=0;j<n2;j++) {
727  O2SCL_IX2(C,i,j)=0.0;
728  }
729  }
730  } else if (beta != 1.0) {
731  for (i=0;i<n1;i++) {
732  for (j=0;j<n2;j++) {
733  O2SCL_IX2(C,i,j)*=beta;
734  }
735  }
736  }
737 
738  if (alpha == 0.0) {
739  return;
740  }
741 
742  TransF=(TransB == o2cblas_ConjTrans) ? o2cblas_Trans : TransB;
743  TransG=(TransA == o2cblas_ConjTrans) ? o2cblas_Trans : TransA;
744 
745  if (TransF == o2cblas_NoTrans && TransG == o2cblas_NoTrans) {
746 
747  /* form C := alpha*A*B+C */
748 
749  for (k=0;k<K;k++) {
750  for (i=0;i<n1;i++) {
751  const double temp=alpha*O2SCL_IX2(B,i,k);
752  if (temp != 0.0) {
753  for (j=0;j<n2;j++) {
754  O2SCL_IX2(C,i,j)+=temp*O2SCL_IX2(A,k,j);
755  }
756  }
757  }
758  }
759 
760  } else if (TransF == o2cblas_NoTrans && TransG == o2cblas_Trans) {
761 
762  /* form C := alpha*A*B'+C */
763 
764  for (i=0;i<n1;i++) {
765  for (j=0;j<n2;j++) {
766  double temp=0.0;
767  for (k=0;k<K;k++) {
768  temp+=O2SCL_IX2(B,i,k)*O2SCL_IX2(A,j,k);
769  }
770  O2SCL_IX2(C,i,j)+=alpha*temp;
771  }
772  }
773 
774  } else if (TransF == o2cblas_Trans && TransG == o2cblas_NoTrans) {
775 
776  for (k=0;k<K;k++) {
777  for (i=0;i<n1;i++) {
778  const double temp=alpha*O2SCL_IX2(B,k,i);
779  if (temp != 0.0) {
780  for (j=0;j<n2;j++) {
781  O2SCL_IX2(C,i,j)+=temp*O2SCL_IX2(A,k,j);
782  }
783  }
784  }
785  }
786 
787  } else if (TransF == o2cblas_Trans && TransG == o2cblas_Trans) {
788 
789  for (i=0;i<n1;i++) {
790  for (j=0;j<n2;j++) {
791  double temp=0.0;
792  for (k=0;k<K;k++) {
793  temp+=O2SCL_IX2(B,k,i)*O2SCL_IX2(A,j,k);
794  }
795  O2SCL_IX2(C,i,j)+=alpha*temp;
796  }
797  }
798 
799  } else {
800  O2SCL_ERR("Unrecognized operation in dgemm().",o2scl::exc_einval);
801  }
802  }
803 
804  return;
805  }
806  //@}
807 
808  /// \name Helper BLAS functions - Subvectors
809  //@{
810  /** \brief Compute \f$ y=\alpha x+y \f$ beginning with index \c ie
811  and ending with index \c N-1
812 
813  This function is used in \ref householder_hv().
814 
815  If \c alpha is identical with zero or <tt>N==ie</tt>, this
816  function will perform no calculations and return without calling
817  the error handler.
818 
819  If <tt>ie</tt> is greater than <tt>N-1</tt> then the error
820  handler will be called if \c O2SCL_NO_RANGE_CHECK is not
821  defined.
822  */
823  template<class vec_t, class vec2_t>
824  void daxpy_subvec(const double alpha, const size_t N, const vec_t &X,
825  vec2_t &Y, const size_t ie) {
826 
827  size_t i;
828 
829  if (alpha == 0.0) return;
830 #if O2SCL_NO_RANGE_CHECK
831 #else
832  if (ie+1>N) {
833  O2SCL_ERR("Invalid index in daxpy_subvec().",o2scl::exc_einval);
834  }
835 #endif
836 
837  const size_t m=(N-ie) % 4;
838 
839  for (i=ie;i<ie+m;i++) {
840  O2SCL_IX(Y,i)+=alpha*O2SCL_IX(X,i);
841  }
842 
843  for (;i+3<N;i+=4) {
844  O2SCL_IX(Y,i)+=alpha*O2SCL_IX(X,i);
845  O2SCL_IX(Y,i+1)+=alpha*O2SCL_IX(X,i+1);
846  O2SCL_IX(Y,i+2)+=alpha*O2SCL_IX(X,i+2);
847  O2SCL_IX(Y,i+3)+=alpha*O2SCL_IX(X,i+3);
848  }
849  }
850 
851  /** \brief Compute \f$ r=x \cdot y \f$ beginning with index \c ie and
852  ending with index \c N-1
853 
854  This function is used in \ref householder_hv().
855 
856  If <tt>ie</tt> is greater than <tt>N-1</tt> then the error
857  handler will be called if \c O2SCL_NO_RANGE_CHECK is not
858  defined.
859  */
860  template<class vec_t, class vec2_t>
861  double ddot_subvec(const size_t N, const vec_t &X, const vec2_t &Y,
862  const size_t ie) {
863  double r=0.0;
864  size_t i;
865 
866 #if O2SCL_NO_RANGE_CHECK
867 #else
868  if (ie+1>N) {
869  O2SCL_ERR("Invalid index in ddot_subvec().",o2scl::exc_einval);
870  }
871 #endif
872 
873  const size_t m=(N-ie) % 4;
874 
875  for (i=ie;i<ie+m;i++) {
876  r+=O2SCL_IX(X,i)*O2SCL_IX(Y,i);
877  }
878 
879  for (;i+3<N;i+=4) {
880  r+=O2SCL_IX(X,i)*O2SCL_IX(Y,i);
881  r+=O2SCL_IX(X,i+1)*O2SCL_IX(Y,i+1);
882  r+=O2SCL_IX(X,i+2)*O2SCL_IX(Y,i+2);
883  r+=O2SCL_IX(X,i+3)*O2SCL_IX(Y,i+3);
884  }
885 
886  return r;
887  }
888 
889  /** \brief Compute the norm of the vector \c X beginning with
890  index \c ie and ending with index \c N-1
891 
892  Used in \ref householder_transform().
893 
894  \note The suffix "2" on the function name indicates that this
895  computes the "2-norm", not that the norm is squared.
896 
897  If <tt>ie</tt> is greater than <tt>N-1</tt> then the error
898  handler will be called if \c O2SCL_NO_RANGE_CHECK is not
899  defined.
900  */
901  template<class vec_t>
902  double dnrm2_subvec(const size_t N, const vec_t &X, const size_t ie) {
903 
904  double scale=0.0;
905  double ssq=1.0;
906 
907 #if O2SCL_NO_RANGE_CHECK
908 #else
909  if (ie+1>N) {
910  O2SCL_ERR("Invalid index in dnrm2_subvec().",o2scl::exc_einval);
911  }
912 #endif
913 
914  if (ie+1==N) {
915  return fabs(O2SCL_IX(X,ie));
916  }
917 
918  for (size_t i=ie;i<N;i++) {
919  const double x=O2SCL_IX(X,i);
920 
921  if (x != 0.0) {
922  const double ax=fabs(x);
923 
924  if (scale<ax) {
925  ssq=1.0+ssq*(scale/ax)*(scale/ax);
926  scale=ax;
927  } else {
928  ssq+=(ax/scale)*(ax/scale);
929  }
930  }
931 
932  }
933 
934  return scale*sqrt(ssq);
935  }
936 
937  /** \brief Compute \f$ x=\alpha x \f$ beginning with index \c ie and
938  ending with index \c N-1
939 
940  This function is used in \ref householder_transform().
941 
942  If <tt>ie</tt> is greater than <tt>N-1</tt> then the error
943  handler will be called if \c O2SCL_NO_RANGE_CHECK is not
944  defined.
945  */
946  template<class vec_t>
947  void dscal_subvec(const double alpha, const size_t N, vec_t &X,
948  const size_t ie) {
949 
950 #if O2SCL_NO_RANGE_CHECK
951 #else
952  if (ie+1>N) {
953  O2SCL_ERR("Invalid index in dscal_subvec().",o2scl::exc_einval);
954  }
955 #endif
956 
957  size_t i;
958  const size_t m=(N-ie) % 4;
959 
960  for (i=ie;i<ie+m;i++) {
961  O2SCL_IX(X,i)*=alpha;
962  }
963 
964  for (;i+3<N;i+=4) {
965  O2SCL_IX(X,i)*=alpha;
966  O2SCL_IX(X,i+1)*=alpha;
967  O2SCL_IX(X,i+2)*=alpha;
968  O2SCL_IX(X,i+3)*=alpha;
969  }
970  }
971  //@}
972 
973  /// \name Helper BLAS functions - Subcolums of a matrix
974  //@{
975  /** \brief Compute \f$ y=\alpha x+y \f$ for a subcolumn of a matrix
976 
977  Given the matrix \c X, define the vector \c x as the column with
978  index \c ic. This function computes \f$ y=\alpha x+y \f$
979  for elements in the vectors \c x and \c y from row \c ir to
980  row \c <tt>M-1</tt> (inclusive). All other elements in \c
981  x and \c y are not referenced.
982 
983  Used in householder_hv_sub().
984  */
985  template<class mat_t, class vec_t>
986  void daxpy_subcol(const double alpha, const size_t M, const mat_t &X,
987  const size_t ir, const size_t ic, vec_t &y) {
988 
989 #if O2SCL_NO_RANGE_CHECK
990 #else
991  if (ir+1>M) {
992  O2SCL_ERR("Invalid index in daxpy_subcol().",o2scl::exc_einval);
993  }
994 #endif
995 
996  if (alpha == 0.0) {
997  return;
998  }
999 
1000  size_t i;
1001  const size_t m=(M-ir) % 4;
1002 
1003  for (i=ir;i<m+ir;i++) {
1004  O2SCL_IX(y,i)+=alpha*O2SCL_IX2(X,i,ic);
1005  }
1006 
1007  for (;i+3<M;i+=4) {
1008  O2SCL_IX(y,i)+=alpha*O2SCL_IX2(X,i,ic);
1009  O2SCL_IX(y,i+1)+=alpha*O2SCL_IX2(X,i+1,ic);
1010  O2SCL_IX(y,i+2)+=alpha*O2SCL_IX2(X,i+2,ic);
1011  O2SCL_IX(y,i+3)+=alpha*O2SCL_IX2(X,i+3,ic);
1012  }
1013 
1014  return;
1015  }
1016 
1017  /** \brief Compute \f$ r=x \cdot y \f$ for a subcolumn of a matrix
1018 
1019  Given the matrix \c X, define the vector \c x as the column with
1020  index \c ic. This function computes \f$ r=x \cdot y \f$
1021  for elements in the vectors \c x and \c y from row \c ir to
1022  row \c <tt>M-1</tt> (inclusive). All other elements in \c
1023  x and \c y are not referenced.
1024 
1025  Used in householder_hv_sub().
1026  */
1027  template<class mat_t, class vec_t>
1028  double ddot_subcol(const size_t M, const mat_t &X, const size_t ir,
1029  const size_t ic, const vec_t &y) {
1030 #if O2SCL_NO_RANGE_CHECK
1031 #else
1032  if (ir+1>M) {
1033  O2SCL_ERR("Invalid index in ddot_subcol().",o2scl::exc_einval);
1034  }
1035 #endif
1036 
1037  double r=0.0;
1038  size_t i;
1039  const size_t m=(M-ir) % 4;
1040 
1041  for (i=ir;i<m+ir;i++) {
1042  r+=O2SCL_IX2(X,i,ic)*O2SCL_IX(y,i);
1043  }
1044 
1045  for (;i+3<M;i+=4) {
1046  r+=O2SCL_IX2(X,i,ic)*O2SCL_IX(y,i);
1047  r+=O2SCL_IX2(X,i+1,ic)*O2SCL_IX(y,i+1);
1048  r+=O2SCL_IX2(X,i+2,ic)*O2SCL_IX(y,i+2);
1049  r+=O2SCL_IX2(X,i+3,ic)*O2SCL_IX(y,i+3);
1050  }
1051 
1052  return r;
1053  }
1054 
1055  /** \brief Compute the norm of a subcolumn of a matrix
1056 
1057  Given the matrix \c A, define the vector \c x as the column with
1058  index \c ic. This function computes the norm of the part of \c x
1059  from row \c ir to row \c <tt>M-1</tt> (inclusive). All other
1060  elements in \c x are not referenced.
1061 
1062  if \c M is zero, then this function silently returns zero
1063  without calling the error handler.
1064 
1065  This function is used in householder_transform_subcol().
1066 
1067  \note The suffix "2" on the function name indicates that
1068  this computes the "2-norm", not that the norm is squared.
1069  */
1070  template<class mat_t>
1071  double dnrm2_subcol(const mat_t &A, const size_t ir, const size_t ic,
1072  const size_t M) {
1073 
1074  double scale=0.0;
1075  double ssq=1.0;
1076  size_t i;
1077 
1078 #if O2SCL_NO_RANGE_CHECK
1079 #else
1080  if (ir+1>M) {
1081  O2SCL_ERR("Invalid index in dnrm2_subcol().",o2scl::exc_einval);
1082  }
1083 #endif
1084 
1085  // Handle the one-element vector case separately
1086  if (ir+1 == M) {
1087  return fabs(O2SCL_IX2(A,ir,ic));
1088  }
1089 
1090  for (i=ir;i<M;i++) {
1091  const double x=O2SCL_IX2(A,i,ic);
1092 
1093  if (x != 0.0) {
1094  const double ax=fabs(x);
1095 
1096  if (scale<ax) {
1097  ssq=1.0+ssq*(scale/ax)*(scale/ax);
1098  scale=ax;
1099  } else {
1100  ssq+=(ax/scale)*(ax/scale);
1101  }
1102  }
1103 
1104  }
1105 
1106  return scale*sqrt(ssq);
1107  }
1108 
1109  /** \brief Compute \f$ x=\alpha x \f$ for a subcolumn of a matrix
1110 
1111  Given the matrix \c A, define the vector \c x as the column with
1112  index \c ic. This function computes \f$ x= \alpha x \f$ for
1113  elements in the vectors \c x from row \c ir to row \c
1114  <tt>M-1</tt> (inclusive). All other elements in \c x are not
1115  referenced.
1116 
1117  Used in householder_transform_subcol().
1118  */
1119  template<class mat_t>
1120  void dscal_subcol(mat_t &A, const size_t ir, const size_t ic,
1121  const size_t M, const double alpha) {
1122 
1123 #if O2SCL_NO_RANGE_CHECK
1124 #else
1125  if (ir+1>M) {
1126  O2SCL_ERR("Invalid index in dscal_subcol().",o2scl::exc_einval);
1127  }
1128 #endif
1129 
1130  size_t i;
1131  const size_t m=(M-ir) % 4;
1132 
1133  for (i=ir;i<m+ir;i++) {
1134  O2SCL_IX2(A,i,ic)*=alpha;
1135  }
1136 
1137  for (;i+3<M;i+=4) {
1138  O2SCL_IX2(A,i,ic)*=alpha;
1139  O2SCL_IX2(A,i+1,ic)*=alpha;
1140  O2SCL_IX2(A,i+2,ic)*=alpha;
1141  O2SCL_IX2(A,i+3,ic)*=alpha;
1142  }
1143 
1144  return;
1145  }
1146 
1147  /** \brief Compute \f$ x=\alpha x \f$ for a subcolumn of a matrix
1148 
1149  Given the matrix \c A, define the vector \c x as the column with
1150  index \c ic. This function computes \f$ x= \alpha x \f$ for
1151  elements in the vectors \c x from row \c ir to row \c
1152  <tt>M-1</tt> (inclusive). All other elements in \c x are not
1153  referenced.
1154 
1155  Used in householder_transform_subcol().
1156  */
1157  template<class mat_t>
1158  double dasum_subcol(mat_t &A, const size_t ir, const size_t ic,
1159  const size_t M) {
1160 
1161 #if O2SCL_NO_RANGE_CHECK
1162 #else
1163  if (ir+1>M) {
1164  O2SCL_ERR("Invalid index in dscal_subcol().",o2scl::exc_einval);
1165  }
1166 #endif
1167 
1168  size_t i;
1169  double r=0.0;
1170  const size_t m=(M-ir) % 4;
1171 
1172  for (i=ir;i<m+ir;i++) {
1173  r+=fabs(O2SCL_IX2(A,i,ic));
1174  }
1175 
1176  for (;i+3<M;i+=4) {
1177  r+=fabs(O2SCL_IX2(A,i,ic));
1178  r+=fabs(O2SCL_IX2(A,i+1,ic));
1179  r+=fabs(O2SCL_IX2(A,i+2,ic));
1180  r+=fabs(O2SCL_IX2(A,i+3,ic));
1181  }
1182 
1183  return r;
1184  }
1185  //@}
1186 
1187  /// \name Helper BLAS functions - Subrows of a matrix
1188  //@{
1189  /** \brief Compute \f$ y=\alpha x+y \f$ for a subrow of a matrix
1190 
1191  Given the matrix \c X, define the vector \c x as the row with
1192  index \c ir. This function computes \f$ y=\alpha x+y \f$ for
1193  elements in the vectors \c x from column \c ic to column \c
1194  <tt>N-1</tt> (inclusive). All other elements in \c x and
1195  \c y are not referenced.
1196 
1197  If <tt>ic</tt> is greater than <tt>N-1</tt> then the error
1198  handler will be called if \c O2SCL_NO_RANGE_CHECK is not
1199  defined.
1200 
1201  Used in householder_hv_sub().
1202  */
1203  template<class mat_t, class vec_t>
1204  void daxpy_subrow(const double alpha, const size_t N, const mat_t &X,
1205  const size_t ir, const size_t ic, vec_t &Y) {
1206 
1207 #if O2SCL_NO_RANGE_CHECK
1208 #else
1209  if (ic+1>N) {
1210  O2SCL_ERR("Invalid index in daxpy_subrow().",o2scl::exc_einval);
1211  }
1212 #endif
1213 
1214  if (alpha == 0.0) {
1215  return;
1216  }
1217 
1218  size_t i;
1219  const size_t m=(N-ic) % 4;
1220 
1221  for (i=ic;i<m+ic;i++) {
1222  O2SCL_IX(Y,i)+=alpha*O2SCL_IX2(X,ir,i);
1223  }
1224 
1225  for (;i+3<N;i+=4) {
1226  O2SCL_IX(Y,i)+=alpha*O2SCL_IX2(X,ir,i);
1227  O2SCL_IX(Y,i+1)+=alpha*O2SCL_IX2(X,ir,i+1);
1228  O2SCL_IX(Y,i+2)+=alpha*O2SCL_IX2(X,ir,i+2);
1229  O2SCL_IX(Y,i+3)+=alpha*O2SCL_IX2(X,ir,i+3);
1230  }
1231 
1232  return;
1233  }
1234 
1235  /** \brief Compute \f$ r=x \cdot y \f$ for a subrow of a matrix
1236 
1237  Given the matrix \c X, define the vector \c x as the row with
1238  index \c ir. This function computes \f$ r=x \cdot y \f$ for
1239  elements in the vectors \c x from column \c ic to column \c
1240  <tt>N-1</tt> (inclusive). All other elements in \c x and
1241  \c y are not referenced.
1242 
1243  If <tt>ic</tt> is greater than <tt>N-1</tt> then the error
1244  handler will be called if \c O2SCL_NO_RANGE_CHECK is not
1245  defined.
1246 
1247  Used in householder_hv_sub().
1248  */
1249  template<class mat_t, class vec_t>
1250  double ddot_subrow(const size_t N, const mat_t &X, const size_t ir,
1251  const size_t ic, const vec_t &Y) {
1252 
1253 #if O2SCL_NO_RANGE_CHECK
1254 #else
1255  if (ic+1>N) {
1256  O2SCL_ERR("Invalid index in ddot_subrow().",o2scl::exc_einval);
1257  }
1258 #endif
1259 
1260  double r=0.0;
1261  size_t i;
1262  const size_t m=(N-ic) % 4;
1263 
1264  for (i=ic;i<m+ic;i++) {
1265  r+=O2SCL_IX2(X,ir,i)*O2SCL_IX(Y,i);
1266  }
1267 
1268  for (;i+3<N;i+=4) {
1269  r+=O2SCL_IX2(X,ir,i)*O2SCL_IX(Y,i);
1270  r+=O2SCL_IX2(X,ir,i+1)*O2SCL_IX(Y,i+1);
1271  r+=O2SCL_IX2(X,ir,i+2)*O2SCL_IX(Y,i+2);
1272  r+=O2SCL_IX2(X,ir,i+3)*O2SCL_IX(Y,i+3);
1273  }
1274 
1275  return r;
1276  }
1277 
1278  /** \brief Compute the norm of a subrow of a matrix
1279 
1280  Given the matrix \c X, define the vector \c x as the row with
1281  index \c ir. This function computes the norm of the part of \c x
1282  from column \c ic to column \c <tt>N-1</tt> (inclusive). All
1283  other elements in \c x are not referenced.
1284 
1285  \note The suffix "2" on the function name indicates that this
1286  computes the "2-norm", not that the norm is squared.
1287  */
1288  template<class mat_t>
1289  double dnrm2_subrow(const mat_t &M, const size_t ir, const size_t ic,
1290  const size_t N) {
1291 
1292  double scale=0.0;
1293  double ssq=1.0;
1294  size_t i;
1295 
1296  if (ic+1==N) {
1297  return fabs(O2SCL_IX2(M,ir,ic));
1298  }
1299 
1300  for (i=ic;i<N;i++) {
1301  const double x=O2SCL_IX2(M,ir,i);
1302 
1303  if (x != 0.0) {
1304  const double ax=fabs(x);
1305 
1306  if (scale<ax) {
1307  ssq=1.0+ssq*(scale/ax)*(scale/ax);
1308  scale=ax;
1309  } else {
1310  ssq+=(ax/scale)*(ax/scale);
1311  }
1312  }
1313 
1314  }
1315 
1316  return scale*sqrt(ssq);
1317  }
1318 
1319  /** \brief Compute \f$ x=\alpha x \f$ for a subrow of a matrix
1320 
1321  Given the matrix \c A, define the vector \c x as the row with
1322  index \c ir. This function computes \f$ x = \alpha x \f$ for
1323  elements in the vectors \c x from column \c ic to column \c
1324  <tt>N-1</tt> (inclusive). All other elements in \c x and
1325  \c y are not referenced.
1326 
1327  If <tt>ic</tt> is greater than <tt>N-1</tt> then the error
1328  handler will be called if \c O2SCL_NO_RANGE_CHECK is not
1329  defined.
1330  */
1331  template<class mat_t>
1332  void dscal_subrow(mat_t &A, const size_t ir, const size_t ic,
1333  const size_t N, const double alpha) {
1334 
1335 #if O2SCL_NO_RANGE_CHECK
1336 #else
1337  if (ic+1>N) {
1338  O2SCL_ERR("Invalid index in dscal_subrow().",o2scl::exc_einval);
1339  }
1340 #endif
1341 
1342  size_t i;
1343  const size_t m=(N-ic) % 4;
1344 
1345  for (i=ic;i<m+ic;i++) {
1346  O2SCL_IX2(A,ir,i)*=alpha;
1347  }
1348 
1349  for (;i+3<N;i+=4) {
1350  O2SCL_IX2(A,ir,i)*=alpha;
1351  O2SCL_IX2(A,ir,i+1)*=alpha;
1352  O2SCL_IX2(A,ir,i+2)*=alpha;
1353  O2SCL_IX2(A,ir,i+3)*=alpha;
1354  }
1355 
1356  return;
1357  }
1358  //@}
1359 
1360 #ifdef DOXYGEN
1361 }
1362 #endif
1363 
void daxpy_subvec(const double alpha, const size_t N, const vec_t &X, vec2_t &Y, const size_t ie)
Compute beginning with index ie and ending with index N-1.
Definition: cblas_base.h:824
double dasum(const size_t N, const vec_t &X)
Compute the absolute sum of vector elements.
Definition: cblas_base.h:111
void dtrsv(const enum o2cblas_order order, const enum o2cblas_uplo Uplo, const enum o2cblas_transpose TransA, const enum o2cblas_diag Diag, const size_t M, const size_t N, const mat_t &A, vec_t &X)
Compute .
Definition: cblas_base.h:335
void daxpy_subcol(const double alpha, const size_t M, const mat_t &X, const size_t ir, const size_t ic, vec_t &y)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:986
void dgemm(const enum o2cblas_order Order, const enum o2cblas_transpose TransA, const enum o2cblas_transpose TransB, const size_t M, const size_t N, const size_t K, const double alpha, const mat_t &A, const mat_t &B, const double beta, mat_t &C)
Compute .
Definition: cblas_base.h:609
bool matrix_is_finite(size_t m, size_t n, mat_t &data)
Test if the first n elements of a matrix are finite.
Definition: cblas_base.h:85
invalid argument supplied by user
Definition: err_hnd.h:59
o2cblas_uplo
Upper- or lower-triangular.
Definition: cblas_base.h:69
double ddot_subrow(const size_t N, const mat_t &X, const size_t ir, const size_t ic, const vec_t &Y)
Compute for a subrow of a matrix.
Definition: cblas_base.h:1250
o2cblas_diag
Unit or generic diagonal.
Definition: cblas_base.h:72
double dnrm2(const size_t N, const vec_t &X)
Compute the norm of the vector X.
Definition: cblas_base.h:183
void dscal_subrow(mat_t &A, const size_t ir, const size_t ic, const size_t N, const double alpha)
Compute for a subrow of a matrix.
Definition: cblas_base.h:1332
void dgemv(const enum o2cblas_order order, const enum o2cblas_transpose TransA, const size_t M, const size_t N, const double alpha, const mat_t &A, const vec_t &X, const double beta, vec2_t &Y)
Compute .
Definition: cblas_base.h:245
Namespace for O2scl CBLAS function templates.
Definition: cblas.h:144
void dtrmv(const enum o2cblas_order Order, const enum o2cblas_uplo Uplo, const enum o2cblas_transpose TransA, const enum o2cblas_diag Diag, const size_t N, const mat_t &A, vec_t &x)
Compute for the triangular matrix A.
Definition: cblas_base.h:477
o2cblas_transpose
Transpose operations.
Definition: cblas_base.h:65
double dasum_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:1158
double ddot(const size_t N, const vec_t &X, const vec2_t &Y)
Compute .
Definition: cblas_base.h:150
double dnrm2_subcol(const mat_t &A, const size_t ir, const size_t ic, const size_t M)
Compute the norm of a subcolumn of a matrix.
Definition: cblas_base.h:1071
double dnrm2_subvec(const size_t N, const vec_t &X, const size_t ie)
Compute the norm of the vector X beginning with index ie and ending with index N-1.
Definition: cblas_base.h:902
double dnrm2_subrow(const mat_t &M, const size_t ir, const size_t ic, const size_t N)
Compute the norm of a subrow of a matrix.
Definition: cblas_base.h:1289
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
double ddot_subvec(const size_t N, const vec_t &X, const vec2_t &Y, const size_t ie)
Compute beginning with index ie and ending with index N-1.
Definition: cblas_base.h:861
void dscal(const double alpha, const size_t N, vec_t &X)
Compute .
Definition: cblas_base.h:217
o2cblas_side
Left or right sided operation.
Definition: cblas_base.h:75
void daxpy(const double alpha, const size_t N, const vec_t &X, vec2_t &Y)
Compute .
Definition: cblas_base.h:125
double ddot_subcol(const size_t M, const mat_t &X, const size_t ir, const size_t ic, const vec_t &y)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:1028
void daxpy_subrow(const double alpha, const size_t N, const mat_t &X, const size_t ir, const size_t ic, vec_t &Y)
Compute for a subrow of a matrix.
Definition: cblas_base.h:1204
o2cblas_order
Matrix order, either column-major or row-major.
Definition: cblas_base.h:62
void dscal_subvec(const double alpha, const size_t N, vec_t &X, const size_t ie)
Compute beginning with index ie and ending with index N-1.
Definition: cblas_base.h:947
void dscal_subcol(mat_t &A, const size_t ir, const size_t ic, const size_t M, const double alpha)
Compute for a subcolumn of a matrix.
Definition: cblas_base.h:1120

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).