Revision 1924
trunk/src/linad99/ad_atlas.cpp (revision 1924)  

1 
/* 

2 
* $Id$ 

3 
* 

4 
* Author: David Fournier 

5 
* Copyright (c) 20082012 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=mmaxmmin+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=mmaxmmin+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; 
Also available in: Unified diff