1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998-2021 Guido Dhondt */
3 /* Copyright (c) 2016 Jean-Marie Verdun */
4
5 /* This program is free software; you can redistribute it and/or */
6 /* modify it under the terms of the GNU General Public License as */
7 /* published by the Free Software Foundation(version 2); */
8 /* */
9
10 /* This program is distributed in the hope that it will be useful, */
11 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
12 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
13 /* GNU General Public License for more details. */
14
15 /* You should have received a copy of the GNU General Public License */
16 /* along with this program; if not, write to the Free Software */
17 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
18
19 #ifdef TAUCS
20
21 #include <stdio.h>
22 #include <math.h>
23 #include <stdlib.h>
24 #include "CalculiX.h"
25 #include "tau.h"
26 #include <taucs.h>
27
28 taucs_ccs_matrix aa[1];
29 void* F=NULL;
30 char* taufactor[]={ "taucs.factor.LLT=true","taucs.factor.mf=true",
31 "taucs.factor.ordering=amd",NULL };
32 char* taufactorooc[]={ "taucs.factor.LLT=true","taucs.factor.mf=true","taucs.factor.ordering=amd","taucs.ooc=true",
33 "taucs.ooc.basename=/tmp/scratch",
34 "taucs.ooc.memory=500000000.0",NULL };
35 char* tausolve[]={ "taucs.factor=false",NULL };
36 char* tausolveooc[]={"taucs.factor=false","taucs.ooc.basename=/tmp/scratch","taucs.ooc.memory=500000000.0",NULL };
37 ITG *irowtau=NULL,*pointtau=NULL;
38 double *autau=NULL;
39 ITG* perm;
40
41
tau_factor(double * ad,double ** aup,double * adb,double * aub,double * sigma,ITG * icol,ITG ** irowp,ITG * neq,ITG * nzs)42 void tau_factor(double *ad, double **aup, double *adb, double *aub,
43 double *sigma,ITG *icol, ITG **irowp,
44 ITG *neq, ITG *nzs){
45
46 ITG i,j,k,l,*irow=NULL;
47 long long ndim;
48 double *au=NULL;
49 double memory_mb = -1.0;
50 ITG mb = -1;
51 ITG ret;
52
53 printf(" Factoring the system of equations using TAUCS\n\n");
54
55 taucs_logfile("stdout");
56
57 au=*aup;
58 irow=*irowp;
59
60 ndim=*neq+*nzs;
61
62 NNEW(autau,double,ndim);
63 NNEW(irowtau,ITG,ndim);
64 NNEW(pointtau,ITG,*neq+1);
65
66 k=ndim;
67 l=*nzs;
68
69 if(*sigma==0.){
70 pointtau[*neq]=ndim;
71 for(i=*neq-1;i>=0;--i){
72 for(j=0;j<icol[i];++j){
73 irowtau[--k]=irow[--l]-1;
74 autau[k]=au[l];
75 }
76 pointtau[i]=--k;
77 irowtau[k]=i;
78 autau[k]=ad[i];
79 }
80 }
81 else{
82 pointtau[*neq]=ndim;
83 for(i=*neq-1;i>=0;--i){
84 for(j=0;j<icol[i];++j){
85 irowtau[--k]=irow[--l]-1;
86 autau[k]=au[l]-*sigma*aub[l];
87 }
88 pointtau[i]=--k;
89 irowtau[k]=i;
90 autau[k]=ad[i]-*sigma*adb[i];
91 }
92 }
93
94 /* defining the matrix */
95
96 aa->n = *neq;
97 aa->m = *neq;
98 aa->flags = TAUCS_SYMMETRIC | TAUCS_LOWER | TAUCS_DOUBLE;
99 aa->colptr = pointtau;
100 aa->rowind = irowtau;
101 aa->values.d = autau;
102
103 if(*neq<5000){
104 taucs_linsolve(aa,&F,0,NULL,NULL,taufactor,NULL);
105 }
106 else{
107 ret = taucs_linsolve(aa,&F,0,NULL,NULL,taufactorooc,NULL);
108
109 printf(" Return Code from Factoring %" ITGFORMAT "\n\n",ret);
110
111 }
112
113 *aup=au;
114 *irowp=irow;
115
116 return;
117 }
118
tau_solve(double * b,ITG * neq)119 void tau_solve(double *b,ITG *neq){
120
121 ITG i;
122 /*static ITG j;*/
123 double *x=NULL;
124 ITG ret;
125
126 NNEW(x,double,*neq);
127
128 if(*neq<5000){
129 taucs_linsolve(aa,&F,1,x,b,tausolve,NULL);
130 }
131 else{
132 ret = taucs_linsolve(aa,&F,1,x,b,tausolveooc,NULL);
133
134 // ret = taucs_ooc_solve_llt(F, x, b);
135
136 printf(" Return Code from Solving %" ITGFORMAT "\n\n",ret);
137
138 // taucs_io_delete(F);
139
140 }
141
142 for(i=0;i<=*neq-1;++i){
143 b[i]=x[i];
144 }
145 SFREE(x);/*
146 if (mb > 0)
147 memory_mb = (double) mb;
148 else
149 memory_mb = ((double) (-mb)) * taucs_available_memory_size()/1048576.0;
150 */
151 /*j++;printf("%" ITGFORMAT "\n",j);*/
152
153 return;
154 }
155
tau_cleanup()156 void tau_cleanup(){
157
158 // taucs_linsolve(NULL,&F,0,NULL,NULL,NULL,NULL);
159 SFREE(pointtau);
160 SFREE(irowtau);
161 SFREE(autau);
162
163 return;
164 }
165
tau(double * ad,double ** aup,double * adb,double * aub,double * sigma,double * b,ITG * icol,ITG ** irowp,ITG * neq,ITG * nzs)166 void tau(double *ad, double **aup, double *adb, double *aub, double *sigma,
167 double *b, ITG *icol, ITG **irowp,
168 ITG *neq, ITG *nzs){
169
170 if(*neq==0) return;
171
172
173 tau_factor(ad,aup,adb,aub,sigma,icol,irowp,
174 neq,nzs);
175
176 tau_solve(b,neq);
177
178 tau_cleanup();
179
180
181 return;
182 }
183
184 #endif
185