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

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