root/src/tau.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. tau_factor
  2. tau_solve
  3. tau_cleanup
  4. tau

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

/* [<][>][^][v][top][bottom][index][help] */