Revision 1924

trunk/src/linad99/ad_atlas.cpp (revision 1924)
1
/*
2
 * $Id$
3
 *
4
 * Author: David Fournier
5
 * Copyright (c) 2008-2012 Regents of the University of California
6
 */
7
/**
8
 * \file
9
 * Description not yet available.
10
 */
11
#include <admodel.h>
12

  
13
#if defined USE_ATLAS
14

  
15
#include <clapack.h>
16

  
17
/**
18
 * Description not yet available.
19
 * \param
20
 */
21
dvector atlas_solve_spd(const dmatrix & M, const dvector & x)
22
{
23
  int mmin=M.indexmin();
24
  int mmax=M.indexmax();
25
  if (mmin != M(mmin).indexmin() ||
26
     mmax != M(mmin).indexmax())
27
  {
28
    cerr << "Matrix not square in dvector atlas_solve_spd"
29
         << endl;
30
    ad_exit(1);
31
  }
32
  if (mmin != x.indexmin() || mmax != x.indexmax())
33
  {
34
    cerr << "Incompatible matrix and vector sizes in dvector atlas_solve_spd"
35
         << endl;
36
    ad_exit(1);
37
  }
38
  dvector v(mmin,mmax);
39
  int sz=mmax-mmin+1;
40
  dvector M1(1,sz*sz);
41
  int ii=1;
42
  int i,j;
43
  for (i=mmin;i<=mmax;i++)
44
  {
45
    for (j=mmin;j<=mmax;j++)
46
    {
47
      M1(ii++)=M(i,j);
48
    }
49
  }
50
  v=x;
51
  double *Ap= &(M1(1));
52
  double *X = &(v(mmin)); const int incX=1;
53

  
54
  const enum CBLAS_ORDER Order=CblasRowMajor;
55
  const enum CBLAS_UPLO Uplo=CblasLower;
56
  const enum CBLAS_TRANSPOSE TransA=CblasNoTrans;
57
  const enum CBLAS_DIAG Diag=CblasNonUnit;
58

  
59
  int retr=clapack_dposv(Order, Uplo, sz,1, Ap, sz, X, sz);
60
  return v;
61
}
62

  
63
/**
64
 * Description not yet available.
65
 * \param
66
 */
67
dvector atlas_solve_spd(const dmatrix & M, const dvector & x, int& ierr)
68
{
69
  int mmin=M.indexmin();
70
  int mmax=M.indexmax();
71
  if (mmin != M(mmin).indexmin() ||
72
     mmax != M(mmin).indexmax())
73
  {
74
    cerr << "Matrix not square in dvector atlas_solve_spd"
75
         << endl;
76
    ad_exit(1);
77
  }
78
  if (mmin != x.indexmin() || mmax != x.indexmax())
79
  {
80
    cerr << "Incompatible matrix and vector sizes in dvector atlas_solve_spd"
81
         << endl;
82
    ad_exit(1);
83
  }
84
  dvector v(mmin,mmax);
85
  int sz=mmax-mmin+1;
86
  dvector M1(1,sz*sz);
87
  int ii=1;
88
  int i,j;
89
  for (i=mmin;i<=mmax;i++)
90
  {
91
    for (j=mmin;j<=mmax;j++)
92
    {
93
      M1(ii++)=M(i,j);
94
    }
95
  }
96
  v=x;
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff