← Back to branch summary

~ubuntu-branches/ubuntu/karmic/hypre/karmic

« back to all changes in this revision

Viewing changes to src/parcsr_block_mv/csr_block_matrix.c

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2009-03-20 11:40:12 UTC
  • mfrom: (4.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20090320114012-132h6ok9w2r6o609
Tags: 2.4.0b-2
Rebuild against new openmpi.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/parcsr_block_mv/csr_block_matrix.c
1
1
/*BHEADER**********************************************************************
2
 
 * Copyright (c) 2006   The Regents of the University of California.
 
2
 * Copyright (c) 2008,  Lawrence Livermore National Security, LLC.
3
3
 * Produced at the Lawrence Livermore National Laboratory.
4
 
 * Written by the HYPRE team. UCRL-CODE-222953.
5
 
 * All rights reserved.
6
 
 *
7
 
 * This file is part of HYPRE (see http://www.llnl.gov/CASC/hypre/).
8
 
 * Please see the COPYRIGHT_and_LICENSE file for the copyright notice, 
9
 
 * disclaimer, contact information and the GNU Lesser General Public License.
10
 
 *
11
 
 * HYPRE is free software; you can redistribute it and/or modify it under the 
12
 
 * terms of the GNU General Public License (as published by the Free Software
13
 
 * Foundation) version 2.1 dated February 1999.
14
 
 *
15
 
 * HYPRE is distributed in the hope that it will be useful, but WITHOUT ANY 
16
 
 * WARRANTY; without even the IMPLIED WARRANTY OF MERCHANTABILITY or FITNESS 
17
 
 * FOR A PARTICULAR PURPOSE.  See the terms and conditions of the GNU General
18
 
 * Public License for more details.
19
 
 *
20
 
 * You should have received a copy of the GNU Lesser General Public License
21
 
 * along with this program; if not, write to the Free Software Foundation,
22
 
 * Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23
 
 *
24
 
 * $Revision: 1.10 $
 
4
 * This file is part of HYPRE.  See file COPYRIGHT for details.
 
5
 *
 
6
 * HYPRE is free software; you can redistribute it and/or modify it under the
 
7
 * terms of the GNU Lesser General Public License (as published by the Free
 
8
 * Software Foundation) version 2.1 dated February 1999.
 
9
 *
 
10
 * $Revision: 1.14 $
25
11
 ***********************************************************************EHEADER*/
26
12
 
27
13
 
28
14
 
 
15
 
29
16
/******************************************************************************
30
17
 *
31
18
 * Member functions for hypre_CSRBlockMatrix class.
374
361
}
375
362
 
376
363
/*--------------------------------------------------------------------------
 
364
 * hypre_CSRBlockMatrixBlockAddAccumulateDiagCheckSign
 
365
 *only add elements of sign*i1 that are negative (sign is size block_size)
 
366
 
 
367
 * (diag(o) = diag(i1) + diag(o)) 
 
368
 *--------------------------------------------------------------------------*/
 
369
int
 
370
hypre_CSRBlockMatrixBlockAddAccumulateDiagCheckSign(double* i1, double* o, int block_size, double *sign)
 
371
{
 
372
   int i;
 
373
   double tmp;
 
374
 
 
375
   for (i = 0; i < block_size; i++)
 
376
   {
 
377
      
 
378
      tmp = i1[i*block_size+i]*sign[i];
 
379
      if (tmp < 0) 
 
380
         o[i*block_size+i] += i1[i*block_size+i];
 
381
   }
 
382
   
 
383
   return 0;
 
384
}
 
385
 
 
386
/*--------------------------------------------------------------------------
 
387
 *  hypre_CSRBlockMatrixComputeSign
 
388
 
 
389
 * o = sign(diag(i1)) 
 
390
 *--------------------------------------------------------------------------*/
 
391
 
 
392
int hypre_CSRBlockMatrixComputeSign(double *i1, double *o, int block_size)
 
393
{
 
394
 
 
395
   int i;
 
396
   
 
397
 
 
398
    for (i = 0; i < block_size; i++)
 
399
   {
 
400
      if (i1[i*block_size+i] < 0)
 
401
         o[i] = -1;
 
402
      else
 
403
         o[i] = 1;
 
404
   }
 
405
 
 
406
   return 0;
 
407
   
 
408
}
 
409
 
 
410
 
 
411
/*--------------------------------------------------------------------------
377
412
 * hypre_CSRBlockMatrixBlockSetScalar
378
413
 * (each entry in block o is set to beta ) 
379
414
 *--------------------------------------------------------------------------*/
451
486
 * hypre_CSRBlockMatrixBlockNorm
452
487
 * (out = norm(data) ) 
453
488
 *
454
 
 *  (note these are not all actually "norms")
 
489
 *  (note: these are not all actually "norms")
455
490
 *  
456
491
 *
457
492
 *--------------------------------------------------------------------------*/
468
503
 
469
504
   switch (norm_type)
470
505
   {
471
 
      
 
506
      case 6: /* sum of all elements in the block  */
 
507
      {
 
508
         for(i = 0; i < sz; i++) 
 
509
         {
 
510
               sum += (data[i]);
 
511
         }
 
512
         break;
 
513
      }
472
514
      case 5: /* one norm  - max col sum*/
473
515
      {
474
516
        
554
596
   
555
597
}
556
598
 
 
599
 
557
600
/*--------------------------------------------------------------------------
558
601
 * hypre_CSRBlockMatrixBlockMultAdd
559
602
 * (o = i1 * i2 + beta * o) 
659
702
}
660
703
 
661
704
/*--------------------------------------------------------------------------
 
705
 * hypre_CSRBlockMatrixBlockMultAddDiagCheckSign
 
706
 * 
 
707
 *  only mult elements if sign*diag(i2) is negative
 
708
 *(diag(o) = diag(i1) * diag(i2) + beta * diag(o)) 
 
709
 *--------------------------------------------------------------------------*/
 
710
int
 
711
hypre_CSRBlockMatrixBlockMultAddDiagCheckSign(double* i1, double* i2, double beta, 
 
712
                                              double* o, int block_size, double *sign)
 
713
{
 
714
   int    i;
 
715
   double tmp;
 
716
   
 
717
 
 
718
   if (beta == 0.0)
 
719
   {
 
720
      for (i = 0; i < block_size; i++)
 
721
      {
 
722
         tmp = i2[i*block_size+i]*sign[i];
 
723
         if (tmp < 0)
 
724
            o[i*block_size + i] = i1[i*block_size + i] * i2[i*block_size + i];
 
725
      }
 
726
   }
 
727
   else if (beta == 1.0)
 
728
   {
 
729
      for(i = 0; i < block_size; i++)
 
730
      {
 
731
         tmp = i2[i*block_size+i]*sign[i];
 
732
         if (tmp < 0)
 
733
            o[i*block_size + i] = o[i*block_size + i] + i1[i*block_size + i] * i2[i*block_size + i];
 
734
      }
 
735
   }
 
736
   else
 
737
   {
 
738
      for(i = 0; i < block_size; i++)
 
739
      {
 
740
         tmp = i2[i*block_size+i]*sign[i];
 
741
         if (tmp < 0)
 
742
            o[i*block_size + i] = beta* o[i*block_size + i] + i1[i*block_size + i] * i2[i*block_size + i];
 
743
      }
 
744
   }
 
745
   return 0;
 
746
}
 
747
 
 
748
 
 
749
 
 
750
/*--------------------------------------------------------------------------
 
751
 * hypre_CSRBlockMatrixBlockMultAddDiag2 (scales cols of il by diag of i2)
 
752
 * ((o) = (i1) * diag(i2) + beta * (o)) 
 
753
 *--------------------------------------------------------------------------*/
 
754
int
 
755
hypre_CSRBlockMatrixBlockMultAddDiag2(double* i1, double* i2, double beta, 
 
756
                                 double* o, int block_size)
 
757
{
 
758
   int    i, j;
 
759
 
 
760
   if (beta == 0.0)
 
761
   {
 
762
 
 
763
      for (i = 0; i < block_size; i++)
 
764
      {
 
765
         for (j = 0; j < block_size; j++)
 
766
         {
 
767
            o[i*block_size + j] =  i1[i*block_size + j] * i2[j*block_size + j];
 
768
            
 
769
         }
 
770
      }
 
771
   }
 
772
   else if (beta == 1.0)
 
773
   {
 
774
      for (i = 0; i < block_size; i++)
 
775
      {
 
776
         for (j = 0; j < block_size; j++)
 
777
         {
 
778
            o[i*block_size + j] =  o[i*block_size + j] +  i1[i*block_size + j] * i2[j*block_size + j];
 
779
            
 
780
         }
 
781
      }
 
782
 
 
783
 
 
784
   }
 
785
   else
 
786
   {
 
787
     for (i = 0; i < block_size; i++)
 
788
      {
 
789
         for (j = 0; j < block_size; j++)
 
790
         {
 
791
            o[i*block_size + j] =  beta * o[i*block_size + j] +  i1[i*block_size + j] * i2[j*block_size + j];
 
792
            
 
793
         }
 
794
      }
 
795
   }
 
796
   return 0;
 
797
}
 
798
/*--------------------------------------------------------------------------
 
799
 * hypre_CSRBlockMatrixBlockMultAddDiag3 (scales cols of il by i2 - 
 
800
                                          whose diag elements are row sums)
 
801
 * ((o) = (i1) * diag(i2) + beta * (o)) 
 
802
 *--------------------------------------------------------------------------*/
 
803
int
 
804
hypre_CSRBlockMatrixBlockMultAddDiag3(double* i1, double* i2, double beta, 
 
805
                                 double* o, int block_size)
 
806
{
 
807
   int    i, j;
 
808
 
 
809
   double *row_sum;
 
810
 
 
811
   row_sum = hypre_CTAlloc(double, block_size);
 
812
   for (i = 0; i < block_size; i++)
 
813
      {
 
814
         for (j = 0; j < block_size; j++)
 
815
         {
 
816
            row_sum[i] +=  i2[i*block_size + j];
 
817
         }
 
818
      }
 
819
   
 
820
   if (beta == 0.0)
 
821
   {
 
822
      for (i = 0; i < block_size; i++)
 
823
      {
 
824
         for (j = 0; j < block_size; j++)
 
825
         {
 
826
            o[i*block_size + j] =  i1[i*block_size + j] * row_sum[j];
 
827
            
 
828
         }
 
829
      }
 
830
   }
 
831
   else if (beta == 1.0)
 
832
   {
 
833
      for (i = 0; i < block_size; i++)
 
834
      {
 
835
         for (j = 0; j < block_size; j++)
 
836
         {
 
837
            o[i*block_size + j] =  o[i*block_size + j] +  i1[i*block_size + j] * row_sum[j];
 
838
            
 
839
         }
 
840
      }
 
841
 
 
842
 
 
843
   }
 
844
   else
 
845
   {
 
846
     for (i = 0; i < block_size; i++)
 
847
      {
 
848
         for (j = 0; j < block_size; j++)
 
849
         {
 
850
            o[i*block_size + j] =  beta * o[i*block_size + j] +  i1[i*block_size + j] * row_sum[j];
 
851
            
 
852
         }
 
853
      }
 
854
   }
 
855
 
 
856
   hypre_TFree(row_sum);
 
857
      
 
858
   return 0;
 
859
}
 
860
/*--------------------------------------------------------------------------
662
861
 * hypre_CSRBlockMatrixBlockMatvec
663
862
 * (ov = alpha* mat * v + beta * ov)
664
863
 * mat is the matrix - size is block_size^2 
1276
1475
   return (ierr);
1277
1476
}
1278
1477
 
 
1478
/*--------------------------------------------------------------------------
 
1479
 * hypre_CSRBlockMatrixBlockInvMultDiag2
 
1480
 * (o = (i1)* diag(i2)^-1) - so this scales the cols of il by 
 
1481
                             the diag entries in i2 
 
1482
 *--------------------------------------------------------------------------*/
 
1483
int
 
1484
hypre_CSRBlockMatrixBlockInvMultDiag2(double* i1, double* i2, double* o, int block_size)
 
1485
{
 
1486
 
 
1487
   int ierr = 0;
 
1488
   int i, j;
 
1489
 
 
1490
   double eps, tmp;
 
1491
   
 
1492
   eps = 1.0e-8;
 
1493
 
 
1494
   for (i = 0; i < block_size; i++)
 
1495
   {
 
1496
      if (fabs(i2[i*block_size + i]) > eps)
 
1497
      {
 
1498
         tmp = 1 / i2[i*block_size + i];
 
1499
      }
 
1500
      else
 
1501
      {
 
1502
         tmp = 1.0;
 
1503
      }
 
1504
      for (j=0; j< block_size; j++)  /* this should be re-written to access by row (not col)! */
 
1505
      {
 
1506
         o[j*block_size + i] = i1[j*block_size + i] *tmp;
 
1507
      }
 
1508
   }
 
1509
   
 
1510
   return (ierr);
 
1511
}
 
1512
 
 
1513
 
 
1514
/*--------------------------------------------------------------------------
 
1515
 * hypre_CSRBlockMatrixBlockInvMultDiag3
 
1516
 * (o = (i1)* diag(i2)^-1) - so this scales the cols of il by 
 
1517
                             the i2 whose diags are row sums 
 
1518
 *--------------------------------------------------------------------------*/
 
1519
int
 
1520
hypre_CSRBlockMatrixBlockInvMultDiag3(double* i1, double* i2, double* o, int block_size)
 
1521
{
 
1522
 
 
1523
   int ierr = 0;
 
1524
   int i, j;
 
1525
   double eps, tmp, row_sum;
 
1526
   
 
1527
   eps = 1.0e-8;
 
1528
 
 
1529
   for (i = 0; i < block_size; i++)
 
1530
   {
 
1531
      /* get row sum of i2, row i */
 
1532
      row_sum = 0.0;
 
1533
      for (j=0; j< block_size; j++)  
 
1534
      {
 
1535
         row_sum += i2[i*block_size + j];   
 
1536
      }
 
1537
      
 
1538
      /* invert */
 
1539
      if (fabs(row_sum) > eps)
 
1540
      {
 
1541
         tmp = 1 / row_sum;
 
1542
      }
 
1543
      else
 
1544
      {
 
1545
         tmp = 1.0;
 
1546
      }
 
1547
      /* scale col of i1 */
 
1548
      for (j=0; j< block_size; j++)  /* this should be re-written to access by row (not col)! */
 
1549
      {
 
1550
         o[j*block_size + i] = i1[j*block_size + i] *tmp;
 
1551
      }
 
1552
   }
 
1553
   
 
1554
   return (ierr);
 
1555
}
 
1556
 
 
1557
 
1279
1558
 
1280
1559
 
1281
1560
/*--------------------------------------------------------------------------