Standard C/C++: fetch the files lib.cpp
and lib.h
. You can make a directory where you store
these files, and eventually its compiled version lib.o. The example here is program1.cpp from chapter 6 and performs the matrix inversion.
// Simple matrix inversion example
#include <iostream>
#include <new>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include "lib.h"
using namespace std;
/* function declarations */
void inverse(double **, int);
void inverse(double **a, int n)
{
int i,j, *indx;
double d, *col, **y;
// allocate space in memory
indx = new int[n];
col = new double[n];
y = (double **) matrix(n, n, sizeof(double));
ludcmp(a, n, indx, &d); // LU decompose a[][]
printf("\n\nLU form of matrix of a[][]:\n");
for(i = 0; i < n; i++) {
printf("\n");
for(j = 0; j < n; j++) {
printf(" a[%2d][%2d] = %12.4E",i, j, a[i][j]);
// find inverse of a[][] by columns
for(j = 0; j < n; j++) {
// initialize right-side of linear equations
for(i = 0; i < n; i++) col[i] = 0.0;
col[j] = 1.0;
lubksb(a, n, indx, col);
// save result in y[][]
for(i = 0; i < n; i++) y[i][j] = col[i];
} //j-loop over columns
// return the inverse matrix in a[][]
for(i = 0; i < n; i++) {
for(j = 0; j < n; j++) a[i][j] = y[i][j];
free_matrix((void **) y); // release local memory
delete [] col;
delete []indx;
} // End: function inverse()
#include <iostream>
#include "armadillo"
using namespace arma;
using namespace std;
int main()
{
mat A = randu<mat>(5,5);
vec b = randu<vec>(5);
A.print("A =");
b.print("b=");
// solve Ax = b
vec x = solve(A,b);
// print x
x.print("x=");
// find LU decomp of A, if needed, P is the permutation matrix
mat L, U;
lu(L,U,A);
// print l
L.print(" L= ");
// print U
U.print(" U= ");
//Check that A = LU
(A-P*L*U).print("Test of LU decomposition");
return 0;
}
Till now we have not paid much attention to speed and possible optimization possibilities inherent in the various compilers. We have compiled and linked as
c++ -c mycode.cpp
c++ -o mycode.exe mycode.o
For Fortran replace with for example gfortran or ifort. This is what we call a flat compiler option and should be used when we develop the code. It produces normally a very large and slow code when translated to machine instructions. We use this option for debugging and for establishing the correct program output because every operation is done precisely as the user specified it.
It is instructive to look up the compiler manual for further instructions by writing
man c++
We have additional compiler options for optimization. These may include procedure inlining where performance may be improved, moving constants inside loops outside the loop, identify potential parallelism, include automatic vectorization or replace a division with a reciprocal and a multiplication if this speeds up the code.
c++ -O3 -c mycode.cpp
c++ -O3 -o mycode.exe mycode.o
This (other options are -O2 or -Ofast) is the recommended option.
It is also useful to profile your program under the development stage. You would then compile with
c++ -pg -O3 -c mycode.cpp
c++ -pg -O3 -o mycode.exe mycode.o
After you have run the code you can obtain the profiling information via
gprof mycode.exe > ProfileOutput
When you have profiled properly your code, you must take out this option as it slows down performance. For memory tests use valgrind. An excellent environment for all these aspects, and much more, is Qt creator.
Adding debugging options is a very useful alternative under the development stage of a program. You would then compile with
c++ -g -O0 -c mycode.cpp
c++ -g -O0 -o mycode.exe mycode.o
This option generates debugging information allowing you to trace for example if an array is properly allocated. Some compilers work best with the no optimization option -O0.
Depending on the compiler, one can add flags which generate code that catches integer overflow errors. The flag -ftrapv does this for the CLANG compiler on OS X operating systems.
In general, irrespective of compiler options, it is useful to
k = n-1;
for (i = 0; i < n; i++){
a[i] = b[i] +c*d;
e = g[k];
}
A better code is
temp = c*d;
for (i = 0; i < n; i++){
a[i] = b[i] + temp;
}
e = g[n-1];
Here we avoid a repeated multiplication inside a loop. Most compilers, depending on compiler flags, identify and optimize such bottlenecks on their own, without requiring any particular action by the programmer. However, it is always useful to single out and avoid code examples like the first one discussed here.
Present CPUs are highly parallel processors with varying levels of parallelism. The typical situation can be described via the following three statements.
Vectorization is a special case of Single Instructions Multiple Data (SIMD) to denote a single instruction stream capable of operating on multiple data elements in parallel. We can think of vectorization as the unrolling of loops accompanied with SIMD instructions.
Vectorization is the process of converting an algorithm that performs scalar operations (typically one operation at the time) to vector operations where a single operation can refer to many simultaneous operations. Consider the following example
for (i = 0; i < n; i++){
a[i] = b[i] + c[i];
}
If the code is not vectorized, the compiler will simply start with the first element and then perform subsequent additions operating on one address in memory at the time.
A SIMD instruction can operate on multiple data elements in one single instruction. It uses the so-called 128-bit SIMD floating-point register. In this sense,vectorization adds some form of parallelism since one instruction is applied to many parts of say a vector.
The number of elements which can be operated on in parallel range from four single-precision floating point data elements in so-called Streaming SIMD Extensions and two double-precision floating-point data elements in Streaming SIMD Extensions 2 to sixteen byte operations in a 128-bit register in Streaming SIMD Extensions 2. Thus, vector-length ranges from 2 to 16, depending on the instruction extensions used and on the data type.
We start with the simple scalar operations given by
for (i = 0; i < n; i++){
a[i] = b[i] + c[i];
}
If the code is not vectorized and we have a 128-bit register to store a 32 bits floating point number, it means that we have \( 3\times 32 \) bits that are not used. For the first element we have
| 0 | 1 | 2 | 3 |
| a[0]= | not used | not used | not used |
| b[0]+ | not used | not used | not used |
| c[0] | not used | not used | not used |
We have thus unused space in our SIMD registers. These registers could hold three additional integers.
If we vectorize the code, we can perform, with a 128-bit register four simultaneous operations, that is we have
for (i = 0; i < n; i+=4){
a[i] = b[i] + c[i];
a[i+1] = b[i+1] + c[i+1];
a[i+2] = b[i+2] + c[i+2];
a[i+3] = b[i+3] + c[i+3];
}
displayed here as
| 0 | 1 | 2 | 3 |
| a[0]= | a[1]= | a[2]= | a[3]= |
| b[0]+ | b[1]+ | b[2]+ | b[3]+ |
| c[0] | c[1] | c[2] | c[3] |
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <iomanip>
#include "time.h"
using namespace std; // note use of namespace
int main (int argc, char* argv[])
{
int i = atoi(argv[1]);
double *a, *b, *c;
a = new double[i];
b = new double[i];
c = new double[i];
for (int j = 0; j < i; j++) {
a[j] = 0.0;
b[j] = cos(j*1.0);
c[j] = sin(j*3.0);
}
clock_t start, finish;
start = clock();
for (int j = 0; j < i; j++) {
a[j] = b[j]+b[j]*c[j];
}
finish = clock();
double timeused = (double) (finish - start)/(CLOCKS_PER_SEC );
cout << setiosflags(ios::showpoint | ios::uppercase);
cout << setprecision(10) << setw(20) << "Time used for vector addition and multiplication=" << timeused << endl;
delete [] a;
delete [] b;
delete [] c;
return 0;
}
We can compile and link without vectorization
c++ -o novec.x vecexample.cpp
vand with vectorization (and additional optimizations)
c++ -O3 -o vec.x vecexample.cpp
The speedup depends on the size of the vectors. In the example here we have run with \( 10^7 \) elements. The example here was run on a PC with ubuntu 14.04 as operating system and an Intel i7-4790 CPU running at 3.60 GHz.
Compphys:~ hjensen$ ./vec.x 10000000
Time used for vector addition = 0.0100000
Compphys:~ hjensen$ ./novec.x 10000000
Time used for vector addition = 0.03000000000
This particular C++ compiler speeds up the above loop operations with a factor of 3. Performing the same operations for \( 10^8 \) elements results only in a factor \( 1.4 \). The result will however vary from compiler to compiler. In general however, with optimization flags like \( -O3 \) or \( -Ofast \), we gain a considerable speedup if our code can be vectorized. Many of these operations can be done automatically by your compiler. These automatic or near automatic compiler techniques improve performance considerably.
Not all loops can be vectorized, as discussed in Intel's guide to vectorization
An important criteria is that the loop counter \( n \) is known at the entry of the loop.
for (int j = 0; j < n; j++) {
a[j] = cos(j*1.0);
}
The variable \( n \) does need to be known at compile time. However, this variable must stay the same for the entire duration of the loop. It implies that an exit statement inside the loop cannot be data dependent.
An exit statement should in general be avoided. If the exit statement contains data-dependent conditions, the loop cannot be vectorized. The following is an example of a non-vectorizable loop
for (int j = 0; j < n; j++) {
a[j] = cos(j*1.0);
if (a[j] < 0 ) break;
}
Avoid loop termination conditions and opt for a single entry loop variable \( n \). The lower and upper bounds have to be kept fixed within the loop.
SIMD instructions perform the same type of operations multiple times. A switch statement leads thus to a non-vectorizable loop since different statemens cannot branch. The following code can however be vectorized since the if statement is implemented as a masked assignment.
for (int j = 0; j < n; j++) {
double x = cos(j*1.0);
if (x > 0 ) {
a[j] = x*sin(j*2.0);
}
else {
a[j] = 0.0;
}
}
These operations can be performed for all data elements but only those elements which the mask evaluates as true are stored. In general, one should avoid branches such as switch, go to, or return statements or if constructs that cannot be treated as masked assignments.
Only the innermost loop of the following example is vectorized
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
a[i][j] += b[i][j];
}
}
The exception is if an original outer loop is transformed into an inner loop as the result of compiler optimizations.
Calls to programmer defined functions ruin vectorization. However, calls to intrinsic functions like \( \sin{x} \), \( \cos{x} \), \( \exp{x} \) etc are allowed since they are normally efficiently vectorized. The following example is fully vectorizable
for (int i = 0; i < n; i++) {
a[i] = log10(i)*cos(i);
}
Similarly, inline functions defined by the programmer, allow for vectorization since the function statements are glued into the actual place where the function is called.
One has to keep in mind that vectorization changes the order of operations inside a loop. A so-called read-after-write statement with an explicit flow dependency cannot be vectorized. The following code
double b = 15.;
for (int i = 1; i < n; i++) {
a[i] = a[i-1] + b;
}
is an example of flow dependency and results in wrong numerical results if vectorized. For a scalar operation, the value \( a[i-1] \) computed during the iteration is loaded into the right-hand side and the results are fine. In vector mode however, with a vector length of four, the values \( a[0] \), \( a[1] \), \( a[2] \) and \( a[3] \) from the previous loop will be loaded into the right-hand side and produce wrong results. That is, we have
a[1] = a[0] + b;
a[2] = a[1] + b;
a[3] = a[2] + b;
a[4] = a[3] + b;
and if the two first iterations are executed at the same by the SIMD instruction, the value of say \( a[1] \) could be used by the second iteration before it has been calculated by the first iteration, leading thereby to wrong results.
On the other hand, a so-called write-after-read statement can be vectorized. The following code
double b = 15.;
for (int i = 1; i < n; i++) {
a[i-1] = a[i] + b;
}
is an example of flow dependency that can be vectorized since no iteration with a higher value of \( i \) can complete before an iteration with a lower value of \( i \). However, such code leads to problems with parallelization.
For C++ programmers it is also worth keeping in mind that an array notation is preferred to the more compact use of pointers to access array elements. The compiler can often not tell if it is safe to vectorize the code.
When dealing with arrays, you should also avoid memory stride, since this slows down considerably vectorization. When you access array element, write for example the inner loop to vectorize using unit stride, that is, access successively the next array element in memory, as shown here
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
a[i][j] += b[i][j];
}
}
We can compile and link without vectorization using the clang c++ (on OSX for example) compiler
clang -o novec.x vecexample.cpp
and with vectorization (and additional optimizations)
clang++ -O3 -Rpass=loop-vectorize -o vec.x vecexample.cpp
The speedup depends on the size of the vectors. In the example here we have run with \( 10^7 \) elements. The example here was run on an IMac17.1 with OSX El Capitan (10.11.4) as operating system and an Intel i5 3.3 GHz CPU.
Compphys:~ hjensen$ ./vec.x 10000000
Time used for norm computation=0.04720500000
Compphys:~ hjensen$ ./novec.x 10000000
Time used for norm computation=0.03311700000
This particular C++ compiler speeds up the above loop operations with a factor of 1.5 Performing the same operations for \( 10^9 \) elements results in a smaller speedup since reading from main memory is required. The non-vectorized code is seemingly faster.
Compphys:~ hjensen$ ./vec.x<PHONE_NUMBER>
Time used for norm computation=58.41391100
Compphys:~ hjensen$ ./novec.x<PHONE_NUMBER>
Time used for norm computation=46.51295300
We will discuss these issues below.
We can compile and link without vectorization with clang compiler
clang++ -o -fno-vectorize novec.x vecexample.cpp
and with vectorization
clang++ -O3 -Rpass=loop-vectorize -o vec.x vecexample.cpp
We can also add vectorization analysis, see for example
clang++ -O3 -Rpass-analysis=loop-vectorize -o vec.x vecexample.cpp
or figure out if vectorization was missed
clang++ -O3 -Rpass-missed=loop-vectorize -o vec.x vecexample.cpp
How do we measure erformance? What is wrong with this code to time a loop?
clock_t start, finish;
start = clock();
for (int j = 0; j < i; j++) {
a[j] = b[j]+b[j]*c[j];
}
finish = clock();
double timeused = (double) (finish - start)/(CLOCKS_PER_SEC );
What happens when the code is executed? The assumption is that the code is ready to execute. But
The first step is to multiply the first row by \( a_0/b_0 \) and subtract it from the second row. This is known as the forward substitution step. We obtain then $$ a_i = 0, $$ $$ b_i = b_i - \frac{a_{i-1}}{b_{i-1}}c_{i-1}, $$ and $$ f_i = f_i - \frac{a_{i-1}}{b_{i-1}}f_{i-1}. $$ At this point the simplified equation, with only an upper triangular matrix takes the form $$ \left( \begin{array}{ccccc} b_0 & c_0 & & & \\ & b_1 & c_1 & & \\ & & \ddots & & \\ & & & b_{m-2} & c_{m-2} \\ & & & & b_{m-1} \end{array} \right)\left( \begin{array}{c} x_0 \\ x_1 \\ \vdots \\ x_{m-2} \\ x_{m-1} \end{array} \right)=\left( \begin{array}{c} f_0 \\ f_1 \\ \vdots \\ f_{m-2} \\ f_{m-1} \\ \end{array} \right) $$
The next step is the backward substitution step. The last row is multiplied by \( c_{N-3}/b_{N-2} \) and subtracted from the second to last row, thus eliminating \( c_{N-3} \) from the last row. The general backward substitution procedure is $$ c_i = 0, $$ and $$ f_{i-1} = f_{i-1} - \frac{c_{i-1}}{b_i}f_i $$ All that remains to be computed is the solution, which is the very straight forward process of $$ x_i = \frac{f_i}{b_i} $$
| Operation | Floating Point |
| Memory Reads | \( 14(N-2) \) |
| Memory Writes | \( 4(N-2) \) |
| Subtractions | \( 3(N-2) \) |
| Multiplications | \( 3(N-2) \) |
| Divisions | \( 4(N-2) \) |
// Forward substitution
// Note that we can simplify by precalculating a[i-1]/b[i-1]
for (int i=1; i < n; i++) {
b[i] = b[i] - (a[i-1]*c[i-1])/b[i-1];
f[i] = g[i] - (a[i-1]*f[i-1])/b[i-1];
}
x[n-1] = f[n-1] / b[n-1];
// Backwards substitution
for (int i = n-2; i >= 0; i--) {
f[i] = f[i] - c[i]*f[i+1]/b[i+1];
x[i] = f[i]/b[i];
}
| Operation | Floating Point |
| Memory Reads | \( 6(N-2) \) |
| Memory Writes | \( 2(N-2) \) |
| Additions | \( 2(N-2) \) |
| Divisions | \( 2(N-2) \) |
// Forward substitution cannot be vectorized
for (int i = 2; i < n; i++) b[i] = b[i] + b[i-1]/d[i-1];
// Backward substitution cannot be vectorized
solution[n-1] = b[n-1]/d[n-1];
for (int i = n-2; i > 0; i--) solution[i] = (b[i]+solution[i+1])/d[i];
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <iomanip>
#include "time.h"
using namespace std; // note use of namespace
int main (int argc, char* argv[])
{
// read in dimension of square matrix
int n = atoi(argv[1]);
double **A, **B;
// Allocate space for the two matrices
A = new double*[n]; B = new double*[n];
for (int i = 0; i < n; i++){
A[i] = new double[n];
B[i] = new double[n];
}
// Set up values for matrix A
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++) {
A[i][j] = cos(i*1.0)*sin(j*3.0);
}
}
clock_t start, finish;
start = clock();
// Then compute the transpose
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++) {
B[i][j]= A[j][i];
}
}
finish = clock();
double timeused = (double) (finish - start)/(CLOCKS_PER_SEC );
cout << setiosflags(ios::showpoint | ios::uppercase);
cout << setprecision(10) << setw(20) << "Time used for setting up transpose of matrix=" << timeused << endl;
// Free up space
for (int i = 0; i < n; i++){
delete[] A[i];
delete[] B[i];
}
delete[] A;
delete[] B;
return 0;
}
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <iomanip>
#include "time.h"
using namespace std; // note use of namespace
int main (int argc, char* argv[])
{
// read in dimension of square matrix
int n = atoi(argv[1]);
double s = 1.0/sqrt( (double) n);
double **A, **B, **C;
// Start timing
clock_t start, finish;
start = clock();
// Allocate space for the two matrices
A = new double*[n]; B = new double*[n]; C = new double*[n];
for (int i = 0; i < n; i++){
A[i] = new double[n];
B[i] = new double[n];
C[i] = new double[n];
}
// Set up values for matrix A and B and zero matrix C
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++) {
double angle = 2.0*M_PI*i*j/ (( double ) n);
A[i][j] = s * ( sin ( angle ) + cos ( angle ) );
B[j][i] = A[i][j];
}
}
// Then perform the matrix-matrix multiplication
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++) {
double sum = 0.0;
for (int k = 0; k < n; k++) {
sum += B[i][k]*A[k][j];
}
C[i][j] = sum;
}
}
// Compute now the Frobenius norm
double Fsum = 0.0;
for (int i = 0; i < n; i++){
for (int j = 0; j < n; j++) {
Fsum += C[i][j]*C[i][j];
}
}
Fsum = sqrt(Fsum);
finish = clock();
double timeused = (double) (finish - start)/(CLOCKS_PER_SEC );
cout << setiosflags(ios::showpoint | ios::uppercase);
cout << setprecision(10) << setw(20) << "Time used for matrix-matrix multiplication=" << timeused << endl;
cout << " Frobenius norm = " << Fsum << endl;
// Free up space
for (int i = 0; i < n; i++){
delete[] A[i];
delete[] B[i];
delete[] C[i];
}
delete[] A;
delete[] B;
delete[] C;
return 0;
}
Why object orientation?
In programming languages, a polymorphic object is an entity, such as a variable or a procedure, that can hold or operate on values of differing types during the program's execution. Because a polymorphic object can operate on a variety of values and types, it can also be used in a variety of programs, sometimes with little or no change by the programmer. The idea of write once, run many, also known as code reusability, is an important characteristic to the programming paradigm known as Object-Oriented Programming (OOP).
OOP describes an approach to programming where a program is viewed as a collection of interacting, but mostly independent software components. These software components are known as objects in OOP and they are typically implemented in a programming language as an entity that encapsulates both data and procedures.
In Fortran a vector or matrix start with \( 1 \), but it is easy
to change a vector so that it starts with zero or even a negative number.
If we have a double precision Fortran vector which starts at \( -10 \) and ends at \( 10 \), we could declare it as
REAL(KIND=8) :: vector(-10:10)
. Similarly, if we want to start at zero and end at 10 we could write
REAL(KIND=8) :: vector(0:10)
.
We have also seen that Fortran allows us to write a matrix addition \( \mathbf{A} = \mathbf{B}+\mathbf{C} \) as
A = B + C
. This means that we have overloaded the addition operator so that it translates this operation into
two loops and an addition of two matrix elements \( a_{ij} = b_{ij}+c_{ij} \).
The way the matrix addition is written is very close to the way we express this relation mathematically. The benefit for the programmer is that our code is easier to read. Furthermore, such a way of coding makes it more likely to spot eventual errors as well.
In Ansi C and C++ arrays start by default from \( i=0 \). Moreover, if we wish to add two matrices we need to explicitely write out the two loops as
for(i=0 ; i < n ; i++) {
for(j=0 ; j < n ; j++) {
a[i][j]=b[i][j]+c[i][j]
However, the strength of C++ is the possibility to define new data types, tailored to some particular problem. Via new data types and overloading of operations such as addition and subtraction, we can easily define sets of operations and data types which allow us to write a matrix addition in exactly the same way as we would do in Fortran. We could also change the way we declare a C++ matrix elements \( a_{ij} \), from \( a[i][j] \) to say \( a(i,j) \), as we would do in Fortran. Similarly, we could also change the default range from \( 0:n-1 \) to \( 1:n \).
To achieve this we need to introduce two important entities in C++ programming, classes and templates.
The function and class declarations are fundamental concepts within C++. Functions are abstractions which encapsulate an algorithm or parts of it and perform specific tasks in a program. We have already met several examples on how to use functions. Classes can be defined as abstractions which encapsulate data and operations on these data. The data can be very complex data structures and the class can contain particular functions which operate on these data. Classes allow therefore for a higher level of abstraction in computing. The elements (or components) of the data type are the class data members, and the procedures are the class member functions.
Classes are user-defined tools used to create multi-purpose software which can be reused by other classes or functions. These user-defined data types contain data (variables) and functions operating on the data.
A simple example is that of a point in two dimensions. The data could be the \( x \) and \( y \) coordinates of a given point. The functions we define could be simple read and write functions or the possibility to compute the distance between two points.
C++ has a class complex in its standard template library (STL). The standard usage in a given function could then look like
// Program to calculate addition and multiplication of two complex numbers
using namespace std;
#include <iostream>
#include <cmath>
#include <complex>
int main()
{
complex<double> x(6.1,8.2), y(0.5,1.3);
// write out x+y
cout << x + y << x*y << endl;
return 0;
where we add and multiply two complex numbers \( x=6.1+\imath 8.2 \) and \( y=0.5+\imath 1.3 \) with the obvious results \( z=x+y=6.6+\imath 9.5 \) and \( z=x\cdot y= -7.61+\imath 12.03 \).
We proceed by splitting our task in three files.
We define first a header file complex.h which contains the declarations of the class. The header file contains the class declaration (data and functions), declaration of stand-alone functions, and all inlined functions, starting as follows
#ifndef Complex_H
#define Complex_H
// various include statements and definitions
#include <iostream> // Standard ANSI-C++ include files
#include <new>
#include ....
class Complex
{...
definition of variables and their character
};
// declarations of various functions used by the class
...
#endif
Next we provide a file complex.cpp where the code and algorithms of
different functions (except inlined functions) declared within the
class are written. The files complex.h
and complex.cpp
are normally
placed in a directory with other classes and libraries we have
defined.
Finally, we discuss here an example of a main program which uses this particular class. An example of a program which uses our complex class is given below. In particular we would like our class to perform tasks like declaring complex variables, writing out the real and imaginary part and performing algebraic operations such as adding or multiplying two complex numbers.
#include "Complex.h"
... other include and declarations
int main ()
{
Complex a(0.1,1.3); // we declare a complex variable a
Complex b(3.0), c(5.0,-2.3); // we declare complex variables b and c
Complex d = b; // we declare a new complex variable d
cout << "d=" << d << ", a=" << a << ", b=" << b << endl;
d = a*c + b/a; // we add, multiply and divide two complex numbers
cout << "Re(d)=" << d.Re() << ", Im(d)=" << d.Im() << endl; // write out of the real and imaginary parts
We include the header file complex.h and define four different complex variables. These are \( a=0.1+\imath 1.3 \), \( b=3.0+\imath 0 \) (note that if you don't define a value for the imaginary part this is set to zero), \( c=5.0-\imath 2.3 \) and \( d=b \). Thereafter we have defined standard algebraic operations and the member functions of the class which allows us to print out the real and imaginary part of a given variable.
class Complex
{
private:
double re, im; // real and imaginary part
public:
Complex (); // Complex c;
Complex (double re, double im = 0.0); // Definition of a complex variable;
Complex (const Complex& c); // Usage: Complex c(a); // equate two complex variables
Complex& operator= (const Complex& c); // c = a; // equate two complex variables, same as previous
....
~Complex () {} // destructor
double Re () const; // double real_part = a.Re();
double Im () const; // double imag_part = a.Im();
double abs () const; // double m = a.abs(); // modulus
friend Complex operator+ (const Complex& a, const Complex& b);
friend Complex operator- (const Complex& a, const Complex& b);
friend Complex operator* (const Complex& a, const Complex& b);
friend Complex operator/ (const Complex& a, const Complex& b);
};
The class is defined via the statement class Complex
. We must first use the key word
class
, which in turn is followed by the user-defined variable name Complex
.
The body of the class, data and functions, is encapsulated within the parentheses {...}
.
Data and specific functions can be private, which means that they cannot be accessed from outside the class.
This means also that access cannot be inherited by other functions outside the class. If we use protected
instead of private
, then data and functions can be inherited outside the class.
The key word public
means that data and functions can be accessed from outside the class.
Here we have defined several functions which can be accessed by functions outside the class.
The declaration friend
means that stand-alone functions can work on privately declared variables of the type
(re, im)
. Data members of a class should be declared as private variables.
The first public function we encounter is a so-called
constructor, which tells how we declare a variable of type Complex
and how this variable is initialized. We have chose three possibilities in the example above:
A declaration like Complex c;
calls the member function Complex()
which can have the following implementation
Complex:: Complex () { re = im = 0.0; }
meaning that it sets the real and imaginary parts to zero. Note the way a member function is defined. The constructor is the first function that is called when an object is instantiated.
Another possibility is
Complex:: Complex () {}
which means that there is no initialization of the real and imaginary parts. The drawback is that a given compiler can then assign random values to a given variable.
A call like Complex a(0.1,1.3);
means that we could call the member function `Complex(double, double)`as
Complex:: Complex (double re_a, double im_a) {
re = re_a; im = im_a; }
The simplest member function are those we defined to extract the real and imaginary part of a variable. Here you have to recall that these are private data, that is they invisible for users of the class. We obtain a copy of these variables by defining the functions
double Complex:: Re () const { return re; }} // getting the real part
double Complex:: Im () const { return im; } // and the imaginary part
Note that we have introduced the declaration const
. What does it mean?
This declaration means that a variabale cannot be changed within a called function.
If we define a variable as
const double p = 3;
and then try to change its value, we will get an error when we
compile our program. This means that constant arguments in functions cannot be changed.
// const arguments (in functions) cannot be changed:
void myfunc (const Complex& c)
{ c.re = 0.2; /* ILLEGAL!! compiler error... */ }
If we declare the function and try to change the value to \( 0.2 \), the compiler will complain by sending an error message.
If we define a function to compute the absolute value of complex variable like
double Complex:: abs () { return sqrt(re*re + im*im);}
without the constant declaration and define thereafter a function
myabs
as
double myabs (const Complex& c)
{ return c.abs(); } // Not ok because c.abs() is not a const func.
the compiler would not allow the c.abs() call in myabs
since Complex::abs
is not a constant member function.
Constant functions cannot change the object's state.
To avoid this we declare the function abs
as
double Complex:: abs () const { return sqrt(re*re + im*im); }
C++ (and Fortran) allow for overloading of operators. That means we
can define algebraic operations on for example vectors or any
arbitrary object. As an example, a vector addition of the type
\( \mathbf{c} = \mathbf{a} + \mathbf{b} \) means that we need to write a small part of
code with a for-loop over the dimension of the array. We would rather
like to write this statement as c = a+b;
as this makes the code much
more readable and close to eventual equations we want to code. To
achieve this we need to extend the definition of operators.
Let us study the declarations in our complex class.
In our main function we have a statement like d = b;
, which means
that we call d.operator= (b)
and we have defined a so-called assignment operator
as a part of the class defined as
Complex& Complex:: operator= (const Complex& c)
{
re = c.re;
im = c.im;
return *this;
}
With this function, statements like
Complex d = b;
or Complex d(b);
make a new object \( d \), which becomes a copy of \( b \).
We can make simple implementations in terms of the assignment
Complex:: Complex (const Complex& c)
{ *this = c; }
which is a pointer to "this object", *this
is the present object,
so *this = c;
means setting the present object equal to \( c \), that is
this->operator= (c);
.
The meaning of the addition operator \( + \) for Complex objects is defined in the
function
Complex operator+ (const Complex& a, const Complex& b); // a+b
The compiler translates c = a + b;
into c = operator+ (a, b);
.
Since this implies the call to function, it brings in an additional overhead. If speed
is crucial and this function call is performed inside a loop, then it is more difficult for a
given compiler to perform optimizations of a loop.
The solution to this is to inline functions. We discussed inlining in chapter 2 of the lecture notes. Inlining means that the function body is copied directly into the calling code, thus avoiding calling the function. Inlining is enabled by the inline keyword
inline Complex operator+ (const Complex& a, const Complex& b)
{ return Complex (a.re + b.re, a.im + b.im); }
Inline functions, with complete bodies must be written in the header file complex.h.
Consider the case c = a + b;
that is, c.operator= (operator+ (a,b));
If operator+
, operator=
and the constructor Complex(r,i)
all
are inline functions, this transforms to
c.re = a.re + b.re;
c.im = a.im + b.im;
by the compiler, i.e., no function calls
The stand-alone function operator+
is a friend of the Complex class
class Complex
{
...
friend Complex operator+ (const Complex& a, const Complex& b);
...
};
so it can read (and manipulate) the private data parts \( re \) and \( im \) via
inline Complex operator+ (const Complex& a, const Complex& b)
{ return Complex (a.re + b.re, a.im + b.im); }
Since we do not need to alter the re and im variables, we can get the values by Re() and Im(), and there is no need to be a friend function
inline Complex operator+ (const Complex& a, const Complex& b)
{ return Complex (a.Re() + b.Re(), a.Im() + b.Im()); }
The multiplication functionality can now be extended to imaginary numbers by the following code
inline Complex operator* (const Complex& a, const Complex& b)
{
return Complex(a.re*b.re - a.im*b.im, a.im*b.re + a.re*b.im);
It will be convenient to inline all functions used by this operator.
To inline the complete expression a*b;
, the constructors and
operator=
must also be inlined. This can be achieved via the following piece of code
inline Complex:: Complex () { re = im = 0.0; }
inline Complex:: Complex (double re_, double im_)
{ ... }
inline Complex:: Complex (const Complex& c)
{ ... }
inline Complex:: operator= (const Complex& c)
{ ... }
// e, c, d are complex
e = c*d;
// first compiler translation:
e.operator= (operator* (c,d));
// result of nested inline functions
// operator=, operator*, Complex(double,double=0):
e.re = c.re*d.re - c.im*d.im;
e.im = c.im*d.re + c.re*d.im;
The definitions operator-
and operator/
follow the same set up.
Finally, if we wish to write to file or another device a complex number using the simple syntax
cout << c;
, we obtain this by defining
the effect of \( < < \) for a Complex object as
ostream& operator<< (ostream& o, const Complex& c)
{ o << "(" << c.Re() << "," << c.Im() << ") "; return o;}
What if we wanted to make a class which takes integers
or floating point numbers with single precision?
A simple way to achieve this is copy and paste our class and replace double
with for
example int
.
C++ allows us to do this automatically via the usage of templates, which
are the C++ constructs for parameterizing parts of
classes. Class templates is a template for producing classes. The declaration consists
of the keyword template
followed by a list of template arguments enclosed in brackets.
We can therefore make a more general class by rewriting our original example as
template<class T>
class Complex
{
private:
T re, im; // real and imaginary part
public:
Complex (); // Complex c;
Complex (T re, T im = 0); // Definition of a complex variable;
Complex (const Complex& c); // Usage: Complex c(a); // equate two complex variables
Complex& operator= (const Complex& c); // c = a; // equate two complex variables, same as previous
We can therefore make a more general class by rewriting our original example as
~Complex () {} // destructor
T Re () const; // T real_part = a.Re();
T Im () const; // T imag_part = a.Im();
T abs () const; // T m = a.abs(); // modulus
friend Complex operator+ (const Complex& a, const Complex& b);
friend Complex operator- (const Complex& a, const Complex& b);
friend Complex operator* (const Complex& a, const Complex& b);
friend Complex operator/ (const Complex& a, const Complex& b);
};
What it says is that Complex
is a parameterized type with \( T \) as a parameter and \( T \) has to be a type such as double
or float.
The class complex is now a class template
and we would define variables in a code as
Complex<double> a(10.0,5.1);
Complex<int> b(1,0);
Member functions of our class are defined by preceding the name of the function with the template
keyword.
Consider the function we defined as Complex:: Complex (double re_a, double im_a)
.
We would rewrite this function as
template<class T>
Complex<T>:: Complex (T re_a, T im_a)
{ re = re_a; im = im_a; }
The member functions are otherwise defined following ordinary member function definitions.
Here follows a very simple first class in the file squared.h
// Not all declarations here
// Class to compute the square of a number
template<class T>
class Squared{
public:
// Default constructor, not used here
Squared(){}
// Overload the function operator()
T operator()(T x){return x*x;}
};
and we would use it as
#include <iostream>
#include "squared.h"
using namespace std;
int main(){
Squared<double> s;
cout << s(3) << endl;
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "/Users/hjensen/Teaching/fys4411/programs/cgm/vectorclass.h"
using namespace std;
Vector ConjugateGradient(Matrix A, Vector b, Vector x0){
int dim = x0.Dimension();
const double tolerance = 1.0e-14;
Vector x(dim),r(dim),v(dim),z(dim);
double c,t,d;
x = x0;
r = b - A*x;
v = r;
c = dot(r,r);
for(int i=0;i<dim;i++){
if(sqrt(dot(v,v))<tolerance){
cerr << "An error has occurred in ConjugateGradient: execution of function terminated" << endl;
break;
}
z = A*v;
t = c/dot(v,z);
x = x + t*v;
r = r - t*z;
d = dot(r,r);
if(sqrt(d) < tolerance)
break;
v = r + (d/c)*v;
c = d;
}
return x;
}
Vector SteepestDescent(Matrix A, Vector b, Vector x0){
int IterMax, i;
int dim = x0.Dimension();
const double tolerance = 1.0e-14;
Vector x(dim),f(dim),z(dim);
double c,alpha,d;
IterMax = 30;
x = x0;
f = A*x-b;
i = 0;
while (i <= IterMax || sqrt(dot(f,f)) < tolerance ){
if(sqrt(dot(f,f))<tolerance){
cerr << "An error has occurred: execution of function terminated" << endl;
break;
}
z = A*f;
c = dot(f,f);
alpha = c/dot(f,z);
x = x - alpha*f;
f = A*x-b;
if(sqrt(dot(f,f)) < tolerance) break;
i++;
}
return x;
}
// Main function begins here
int main(int argc, char * argv[]){
int dim = 2;
Vector x(dim),xsd(dim), b(dim),x0(dim);
Matrix A(dim,dim);
// Set our initial guess
x0(0) = x0(1) = 0;
// Set the matrix
A(0,0) = 3; A(1,0) = 2; A(0,1) = 2; A(1,1) = 6;
cout << "The Matrix A that we are using: " << endl;
A.Print();
cout << endl;
Vector y(dim);
y(0) = 2.;
y(1) = -2.;
cout << "The exact solution is: " << endl;
y.Print();
cout << endl;
b = A*y;
cout << "The right hand side, b, of the expression Ax=b: " << endl;
b.Print();
cout << endl;
x = ConjugateGradient(A,b,x0);
xsd = SteepestDescent(A,b,x0);
cout << "The approximate solution using Conjugate Gradient is: " << endl;
x.Print();
cout << endl;
cout << "The approximate solution using Steepest Descent is: " << endl;
xsd.Print();
cout << endl;
}
#ifndef _vectorclass
#define _vectorclass
#include <cmath>
#include <iostream>
using namespace std;
class Point;
class Vector;
class Matrix;
/********************************/
/* Point Class */
/********************************/
class Point{
private:
int dimension;
double *data;
public:
Point(int dim);
Point(const Point& v);
~Point();
int Dimension() const;
//************************
// User Defined Operators
//************************
int operator==(const Point& v) const;
int operator!=(const Point& v) const;
Point & operator=(const Point& v);
double operator()(const int i) const;
double& operator()(const int i);
void Print() const;
};
/********************************/
/* Vector Class */
/********************************/
class Vector{
private:
int dimension;
double *data;
public:
Vector();
Vector(int dim);
Vector(const Vector& v);
Vector(int col, const Matrix &A);
~Vector();
void Initialize(int dim);
int Dimension() const;
double Length(); /* Euclidean Norm of the Vector */
void Normalize();
double Norm_l1();
double Norm_l2();
double Norm_linf();
double MaxMod();
double ElementofMaxMod();
int MaxModindex();
//************************
// User Defined Operators
//************************
int operator==(const Vector& v) const;
int operator!=(const Vector& v) const;
Vector & operator=(const Vector& v);
double operator()(const int i) const;
double& operator()(const int i);
void Print() const;
void Initialize(double a);
void Initialize(double *v);
};
/********************************/
/* Matrix Class */
/********************************/
class Matrix {
private:
int rows, columns;
double **data;
public:
Matrix(int dim);
Matrix(int rows1, int columns1);
Matrix(const Matrix& m);
Matrix(int num_vectors, const Vector * q);
Matrix(int rows1, int columns1, double **rowptrs);
~Matrix();
int Rows() const;
int Columns() const;
double ** GetPointer();
void GetColumn(int col, Vector &x);
void GetColumn(int col, Vector &x, int rowoffset);
void PutColumn(int col, const Vector &x);
double Norm_l1();
double Norm_linf();
//************************
// User Defined Operators
//************************
Matrix& operator=(const Matrix& m);
double operator()(const int i, const int j) const;
double& operator()(const int i, const int j);
double MaxModInRow(int row);
double MaxModInRow(int row, int starting_column);
int MaxModInRowindex(int row);
int MaxModInRowindex(int row, int starting_column);
double MaxModInColumn(int column);
double MaxModInColumn(int column, int starting_row);
int MaxModInColumnindex(int column);
int MaxModInColumnindex(int column, int starting_row);
void RowSwap(int row1, int row2);
void Print() const;
};
/********************************/
/* Operator Declarations */
/********************************/
// Unitary operator -
Vector operator-(const Vector& v);
// Binary operator +,-
Vector operator+(const Vector& v1, const Vector& v2);
Vector operator-(const Vector& v1, const Vector& v2);
// Vector Scaling (multiplication by a scaler : defined commutatively)
Vector operator*(const double s, const Vector& v);
Vector operator*(const Vector& v, const double s);
// Vector Scaling (division by a scaler)
Vector operator/(const Vector& v, const double s);
Vector operator*(const Matrix& A, const Vector& x);
/********************************/
/* Function Declarations */
/********************************/
int min_dimension(const Vector& u, const Vector& v);
double dot(const Vector& u, const Vector& v);
double dot(int N, double *a, double *b);
double dot(int N, const Vector &u, const Vector &v);
void Swap(double &a, double &b);
double Sign(double x);
/* Misc. useful functions to have */
double log2(double x);
double GammaF(double x);
int Factorial(int n);
double ** CreateMatrix(int m, int n);
void DestroyMatrix(double ** mat, int m, int n);
int ** ICreateMatrix(int m, int n);
void IDestroyMatrix(int ** mat, int m, int n);
#endif
#include "vectorclass.h"
Point::Point(int dim){
dimension = dim;
data = new double[dimension];
for(int i=0;i<dimension;i++)
data[i] = 0.0;
}
Point::Point(const Point &v){
dimension = v.Dimension();
data = new double[dimension];
for(int i=0;i<dimension;i++)
data[i] = v.data[i];
}
Point::~Point(){
dimension = 0;
delete[] data;
data = NULL;
}
int Point::Dimension() const{
return(dimension);
}
double Point::operator()(const int i) const{
if(i>=0 && i<dimension)
return data[i];
cerr << "Point::Invalid index " << i << " for Point of dimension " << dimension << endl;
return(0);
}
double& Point::operator()(const int i){
if(i>=0 && i<dimension)
return data[i];
cerr << "Point::Invalid index " << i << " for Point of dimension " << dimension << endl;
return(data[0]);
}
Point& Point::operator=(const Point &v) {
dimension = v.Dimension();
for(int i=0;i<dimension;i++)
data[i] = v.data[i];
return *this;
};
void Point::Print() const{
cout << endl;
cout << "[ ";
if(dimension>0)
cout << data[0];
for(int i=1;i<dimension;i++)
cout << "; " << data[i];
cout << " ]" << endl;
}
Vector::Vector(){
dimension = 0;
data = NULL;
}
Vector::Vector(int dim){
dimension = dim;
data = new double[dimension];
for(int i=0;i<dimension;i++)
data[i] = 0.0;
}
Vector::Vector(const Vector &v){
dimension = v.Dimension();
data = new double[dimension];
for(int i=0;i<dimension;i++)
data[i] = v.data[i];
}
Vector::Vector(int col, const Matrix &A){
dimension = A.Rows();
data = new double[dimension];
for(int i=0;i<A.Rows();i++)
data[i] = A(i,col);
}
Vector::~Vector(){
dimension = 0;
delete[] data;
data = NULL;
}
void Vector::Initialize(int dim){
if(dimension!=0)
delete[] data;
dimension = dim;
data = new double[dimension];
for(int i=0;i<dimension;i++)
data[i] = 0.0;
}
int Vector::Dimension() const{
return(dimension);
}
double Vector::operator()(const int i) const{
if(i>=0 && i<dimension)
return data[i];
cerr << "Vector::Invalid index " << i << " for Vector of dimension " << dimension << endl;
return(0);
}
double& Vector::operator()(const int i){
if(i>=0 && i<dimension)
return data[i];
cerr << "Vector::Invalid index " << i << " for Vector of dimension " << dimension << endl;
return(data[0]);
}
Vector& Vector::operator=(const Vector &v) {
dimension = v.Dimension();
for(int i=0;i<dimension;i++)
data[i] = v.data[i];
return *this;
};
void Vector::Print() const{
cout << endl;
cout << "[ ";
if(dimension>0)
cout << data[0];
for(int i=1;i<dimension;i++)
cout << "; " << data[i];
cout << " ]" << endl;
}
double Vector::Norm_l1(){
double sum = 0.0;
for(int i=0;i<dimension;i++)
sum += fabs(data[i]);
return(sum);
}
double Vector::Norm_l2(){
double sum = 0.0;
for(int i=0;i<dimension;i++)
sum += data[i]*data[i];
return(sqrt(sum));
}
void Vector::Normalize(){
double tmp = 1.0/Norm_l2();
for(int i=0;i<dimension;i++)
data[i] = data[i]*tmp;
}
double Vector::Norm_linf(){
double maxval = 0.0,tmp;
for(int i=0;i<dimension;i++){
tmp = fabs(data[i]);
maxval = (maxval > tmp)?maxval:tmp;
}
return(maxval);
}
double Vector::MaxMod(){
double maxm = -1.0e+10;
for(int i=0; i<dimension; i++)
maxm = (maxm > fabs(data[i]))?maxm:fabs(data[i]);
return maxm;
}
double Vector::ElementofMaxMod(){
return(data[MaxModindex()]);
}
int Vector::MaxModindex(){
double maxm = -1.0e+10;
int maxmindex = 0;
for(int i=0; i<dimension; i++){
if(maxm<fabs(data[i])){
maxm = fabs(data[i]);
maxmindex = i;
}
}
return maxmindex;
}
void Vector::Initialize(double a){
for(int i=0; i<dimension; i++)
data[i] = a;
}
void Vector::Initialize(double *v){
for(int i=0; i<dimension; i++)
data[i] = v[i];
}
Matrix::Matrix(int dim){
rows = dim;
columns = dim;
data = new double*[rows];
for(int i=0;i<rows;i++){
data[i] = new double[columns];
for(int j=0;j<columns;j++)
data[i][j] = 0.0;
}
}
Matrix::Matrix(int rows1, int columns1){
rows = rows1;
columns = columns1;
data = new double*[rows];
for(int i=0;i<rows;i++){
data[i] = new double[columns];
for(int j=0;j<columns;j++)
data[i][j] = 0.0;
}
}
Matrix::Matrix(const Matrix& m){
rows = m.rows;
columns = m.columns;
data = new double*[rows];
for(int i=0;i<rows;i++){
data[i] = new double[columns];
for(int j=0; j<columns; j++)
data[i][j] = m.data[i][j];
}
}
Matrix::Matrix(int num_Vectors, const Vector * q){
rows = q[0].Dimension();
columns = num_Vectors;
data = new double*[rows];
for(int i=0;i<rows;i++){
data[i] = new double[columns];
for(int j=0; j<columns; j++)
data[i][j] = q[j](i);
}
}
Matrix::Matrix(int rows1, int columns1, double **rowptrs){
rows = rows1;
columns = columns1;
data = new double*[rows];
for(int i=0;i<rows;i++)
data[i] = rowptrs[i];
}
Matrix::~Matrix(){
for(int i=0;i<rows;i++)
delete[] data[i];
rows = 0;
columns = 0;
delete[] data;
}
int Matrix::Rows() const{
return(rows);
}
int Matrix::Columns() const{
return(columns);
}
double **Matrix::GetPointer(){
return(data);
}
void Matrix::GetColumn(int col, Vector &x){
x.Initialize(0.0);
for(int i=0;i<rows;i++)
x(i) = data[i][col];
}
void Matrix::GetColumn(int col, Vector &x, int rowoffset){
x.Initialize(0.0);
for(int i=0;i<rows-rowoffset;i++)
x(i) = data[i+rowoffset][col];
}
void Matrix::PutColumn(int col, const Vector &x){
for(int i=0;i<rows;i++)
data[i][col] = x(i);
}
double Matrix::Norm_linf(){
double maxval = 0.0,sum;
for(int i=0;i<rows;i++){
sum = 0.0;
for(int j=0;j<columns;j++)
sum += fabs(data[i][j]);
maxval = (maxval > sum)?maxval:sum;
}
return(maxval);
}
double Matrix::Norm_l1(){
double maxval = 0.0,sum;
for(int j=0;j<columns;j++){
sum = 0.0;
for(int i=0;i<rows;i++)
sum += fabs(data[i][j]);
maxval = (maxval > sum)?maxval:sum;
}
return(maxval);
}
Matrix& Matrix::operator=(const Matrix &m){
if( (rows == m.rows) && (columns == m.columns)){
for(int i=0; i<rows; i++)
for(int j=0;j<columns;j++){
data[i][j] = m.data[i][j];
}
}
else
cerr << "Matrix Error: Cannot equate matrices of different sizes\n";
return *this;
}
double Matrix::operator()(const int i, const int j) const {
if( (i>=0) && (j>=0) && (i<rows) && (j<columns))
return(data[i][j]);
else
cerr << "Matrix Error: Invalid Matrix indices (" << i << "," << j <<
"), for Matrix of size " << rows << " X " << columns << endl;
return((double)0);
}
double& Matrix::operator()(const int i, const int j) {
if( (i>=0) && (j>=0) && (i<rows) && (j<columns))
return(data[i][j]);
else
cerr << "Matrix Error: Invalid Matrix indices (" << i << "," << j <<
"), for Matrix of size " << rows << " X " << columns << endl;;
return(data[0][0]);
}
void Matrix::Print() const{
cout << endl;
cout << "[ ";
for(int i=0;i<rows;i++){
cout << data[i][0];
for(int j=1;j<columns;j++)
cout << " " << data[i][j];
if(i!=(rows-1))
cout << ";\n";
}
cout << " ]" << endl;
}
double Matrix::MaxModInRow(int row){
double maxv = -1.0e+10;
for(int i=0;i<columns;i++)
maxv = (fabs(data[row][i])>maxv)?fabs(data[row][i]):maxv;
return maxv;
}
double Matrix::MaxModInRow(int row, int starting_column){
double maxv = -1.0e+10;
for(int i=starting_column;i<columns;i++)
maxv = (fabs(data[row][i])>maxv)?fabs(data[row][i]):maxv;
return maxv;
}
int Matrix::MaxModInRowindex(int row){
int maxvindex = 0;
double maxv = -1.0e+10;
for(int i=0;i<columns;i++){
if(maxv < fabs(data[row][i])){
maxv = fabs(data[row][i]);
maxvindex = i;
}
}
return maxvindex;
}
int Matrix::MaxModInRowindex(int row, int starting_column){
int maxvindex = 0;
double maxv = -1.0e+10;
for(int i=starting_column;i<columns;i++){
if(maxv < fabs(data[row][i])){
maxv = fabs(data[row][i]);
maxvindex = i;
}
}
return maxvindex;
}
double Matrix::MaxModInColumn(int column){
double maxv = -1.0e+10;
for(int i=0;i<rows;i++)
maxv = (fabs(data[i][column])>maxv)?fabs(data[i][column]):maxv;
return maxv;
}
double Matrix::MaxModInColumn(int column, int starting_row){
double maxv = -1.0e+10;
for(int i=starting_row;i<rows;i++)
maxv = (fabs(data[i][column])>maxv)?fabs(data[i][column]):maxv;
return maxv;
}
int Matrix::MaxModInColumnindex(int column){
int maxvindex = 0;
double maxv = -1.0e+10;
for(int i=0;i<rows;i++){
if(maxv < fabs(data[i][column])){
maxv = fabs(data[i][column]);
maxvindex = i;
}
}
return maxvindex;
}
int Matrix::MaxModInColumnindex(int column, int starting_column){
int maxvindex = 0;
double maxv = -1.0e+10;
for(int i=starting_column;i<rows;i++){
if(maxv < fabs(data[i][column])){
maxv = fabs(data[i][column]);
maxvindex = i;
}
}
return maxvindex;
}
void Matrix::RowSwap(int row1, int row2){
double * tmp = data[row1];
data[row1] = data[row2];
data[row2] = tmp;
}
/****************************************************************/
/* Operator Definitions */
/****************************************************************/
Vector operator-(const Vector& v){
Vector x(v.Dimension());
for(int i=0;i<v.Dimension();i++)
x(i) = -v(i);
return x;
}
Vector operator+(const Vector& v1, const Vector& v2){
int min_dim = min_dimension(v1,v2);
Vector x(min_dim);
for(int i=0;i<min_dim;i++)
x(i) = v1(i) + v2(i);
return x;
}
Vector operator-(const Vector& v1, const Vector& v2){
int min_dim = min_dimension(v1,v2);
Vector x(min_dim);
for(int i=0;i<min_dim;i++)
x(i) = v1(i) - v2(i);
return x;
}
Vector operator/(const Vector& v, const double s) {
Vector x(v.Dimension());
for(int i=0;i<v.Dimension();i++)
x(i) = v(i)/s;
return x;
}
Vector operator*(const double s, const Vector &v) {
Vector x(v.Dimension());
for(int i=0;i<v.Dimension();i++)
x(i) = s*v(i);
return x;
}
Vector operator*(const Vector& v, const double s) {
Vector x(v.Dimension());
for(int i=0;i<v.Dimension();i++)
x(i) = s*v(i);
return x;
}
Vector operator*(const Matrix& A, const Vector& x){
int rows = A.Rows(), columns = A.Columns();
int dim = x.Dimension();
Vector b(dim);
if(columns != dim){
cerr << "Invalid dimensions given in matrix-vector multiply" << endl;
return(b);
}
for(int i=0;i<rows;i++){
b(i) = 0.0;
for(int j=0;j<columns;j++){
b(i) += A(i,j)*x(j);
}
}
return b;
}
/****************************************************************/
/* Function Definitions */
/****************************************************************/
int min_dimension(const Vector& v1, const Vector& v2){
int min_dim = (v1.Dimension()<v2.Dimension())?v1.Dimension():v2.Dimension();
return(min_dim);
}
double dot(const Vector& u, const Vector& v){
double sum = 0.0;
int min_dim = min_dimension(u,v);
for(int i=0;i<min_dim;i++)
sum += u(i)*v(i);
return sum;
}
double dot(int N, const Vector& u, const Vector& v){
double sum = 0.0;
for(int i=0;i<N;i++)
sum += u(i)*v(i);
return sum;
}
double dot(int N, double *a, double *b){
double sum = 0.0;
for(int i=0;i<N;i++)
sum += a[i]*b[i];
return sum;
}
/*******************************/
/* Log base 2 of a number */
/*******************************/
double log2(double x){
return(log(x)/log(2.0));
}
void Swap(double &a, double &b){
double tmp = a;
a = b;
b = tmp;
}
double Sign(double x){
double xs;
xs = (x>=0.0)?1.0:-1.0;
return xs;
}
//GammaF function valid for x integer, or x (integer+0.5)
double GammaF(double x){
double gamma = 1.0;
if (x == -0.5)
gamma = -2.0*sqrt(M_PI);
else if (!x) return gamma;
else if ((x-(int)x) == 0.5){
int n = (int) x;
double tmp = x;
gamma = sqrt(M_PI);
while(n--){
tmp -= 1.0;
gamma *= tmp;
}
}
else if ((x-(int)x) == 0.0){
int n = (int) x;
double tmp = x;
while(--n){
tmp -= 1.0;
gamma *= tmp;
}
}
return gamma;
}
int Factorial(int n){
int value=1;
for(int i=n;i>0;i--)
value = value*i;
return value;
}
double ** CreateMatrix(int m, int n){
double ** mat;
mat = new double*[m];
for(int i=0;i<m;i++){
mat[i] = new double[n];
for(int j=0;j<m;j++)
mat[i][j] = 0.0;
}
return mat;
}
int ** ICreateMatrix(int m, int n){
int ** mat;
mat = new int*[m];
for(int i=0;i<m;i++){
mat[i] = new int[n];
for(int j=0;j<m;j++)
mat[i][j] = 0;
}
return mat;
}
void DestroyMatrix(double ** mat, int m, int n){
for(int i=0;i<m;i++)
delete[] mat[i];
delete[] mat;
}
void IDestroyMatrix(int ** mat, int m, int n){
for(int i=0;i<m;i++)
delete[] mat[i];
delete[] mat;
}


This is our handsome Snoopy son, we're very proud of him. Thank you Snoopy and Dusty for such a beautiful boy. Dustin is so much like Snoopy. Several of Snoopy's American friends and fans have inquired about him. Dustin's dam Dusty) is a daughter of Am Can Ch Boradors Exch at Devonsleigh, who was actually a littermate to Am Ch Borador By George (Gordy). Dusty's mother was from our old, proven sound Lindenhall/Staghorn dogs. Dustin is a solid boy with lots of leg under him, good turn of stifle and orthopedic soundness, really nice front angles to match the rear and nice forechest, good genetic clearances and most wonderful of all, he has his sire's wonderful temperament. Dustin is available at stud to approved females. With his pedigree and clearances and his own excellent health and soundness, we believe he has a lot to offer.


{% extends "global/Page.html" %} {% load staticfiles otree_tags %} {% block title %}
These are the points you earned in each part. {{ player.score }} and {{ player.rank }} and {{ player.given_type }}


sparklyr provides bindings to Spark’s distributed machine learning library. In particular, sparklyr allows you to access the machine learning routines provided by the spark.ml package. Together with sparklyr’s dplyr interface, you can easily create and tune machine learning workflows on Spark, orchestrated entirely within R.
sparklyr provides three families of functions that you can use with Spark machine learning:
ml_*
)ft_*
)sdf_*
)An analytic workflow with sparklyr might be composed of the following stages. For an example see Example Workflow.
sdf_*
and ft_*
family of functions to generate new columns, or partition your data set,ml_*
family of functions to model your data,Spark’s machine learning library can be accessed from sparklyr through the ml_*
set of functions:
| Function | Description |
|---|---|
ml_kmeans |
K-Means Clustering |
ml_linear_regression |
Linear Regression |
ml_logistic_regression |
Logistic Regression |
ml_survival_regression |
Survival Regression |
ml_generalized_linear_regression |
Generalized Linear Regression |
ml_decision_tree |
Decision Trees |
ml_random_forest |
Random Forests |
ml_gradient_boosted_trees |
Gradient-Boosted Trees |
ml_pca |
Principal Components Analysis |
ml_naive_bayes |
Naive-Bayes |
ml_multilayer_perceptron |
Multilayer Perceptron |
ml_lda |
Latent Dirichlet Allocation |
ml_one_vs_rest |
One vs Rest |
The ml_*
functions take the arguments response
and features
. But features
can also be a formula with main effects (it currently does not accept interaction terms). The intercept term can be omitted by using -1
.
# Equivalent statements
ml_linear_regression(z ~ -1 + x + y)
ml_linear_regression(intercept = FALSE, response = "z", features = c("x", "y"))
A model is often fit not on a dataset as-is, but instead on some transformation of that dataset. Spark provides feature transformers, facilitating many common transformations of data within a Spark DataFrame, and sparklyr exposes these within the ft_*
family of functions. These routines generally take one or more input columns, and generate a new output column formed as a transformation of those columns.
| Function | Description |
|---|---|
ft_binarizer |
Threshold numerical features to binary (0/1) feature |
ft_bucketizer |
Bucketizer transforms a column of continuous features to a column of feature buckets |
ft_discrete_cosine_transform |
Transforms a length NN real-valued sequence in the time domain into another length NN real-valued sequence in the frequency domain |
ft_elementwise_product |
Multiplies each input vector by a provided weight vector, using element-wise multiplication. |
ft_index_to_string |
Maps a column of label indices back to a column containing the original labels as strings |
ft_quantile_discretizer |
Takes a column with continuous features and outputs a column with binned categorical features |
ft_sql_transformer |
Implements the transformations which are defined by a SQL statement |
ft_string_indexer |
Encodes a string column of labels to a column of label indices |
ft_vector_assembler |
Combines a given list of columns into a single vector column |
Functions for interacting with Spark ML model fits.
| Function | Description |
|---|---|
ml_binary_classification_eval |
Calculates the area under the curve for a binary classification model. |
ml_classification_eval |
Calculates performance metrics (i.e. f1, precision, recall, weightedPrecision, weightedRecall, and accuracy) for binary and multiclass classification model. |
ml_tree_feature_importance |
Calculates variable importance for decision trees (i.e. decision trees, random forests, gradient boosted trees). |
ml_saveload |
Save and load model fits. For use with scoring models across platforms (e.g. using a model as an estimator in a Spark application). These functions are currently experimental and not yet ready for production use. |
Functions for creating custom wrappers to other Spark ML algorithms.
| Function | Description |
|---|---|
ensure |
Enforces Specific Structure for R Objects. |
ml_create_dummy_variables |
Given a column in a Spark DataFrame, generate a new Spark DataFrame containing dummy variable columns. |
ml_model |
Creates an ML Model Object. |
ml_options |
Provides Options for Spark.ML Routines. |
ml_prepare_dataframe |
Prepares a Spark DataFrame for Spark ML Routines. |
ml_prepare_response_features_intercept |
Pre-process / normalize the inputs typically passed to a Spark ML routine. |
We will use the iris
data set to examine a handful of learning algorithms and transformers. The iris data set measures attributes for 150 flowers in 3 different species of iris.
library(sparklyr)
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
sc <- spark_connect(master = "local")
## * Using Spark: 2.1.0
iris_tbl <- copy_to(sc, iris, "iris", overwrite = TRUE)
iris_tbl
## # Source: table<iris> [?? x 5]
## # Database: spark_connection
## Sepal_Length Sepal_Width Petal_Length Petal_Width Species
## <dbl> <dbl> <dbl> <dbl> <chr>
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
## 7 4.6 3.4 1.4 0.3 setosa
## 8 5.0 3.4 1.5 0.2 setosa
## 9 4.4 2.9 1.4 0.2 setosa
## 10 4.9 3.1 1.5 0.1 setosa
## # ... with 140 more rows
Use Spark’s K-means clustering to partition a dataset into groups. K-means clustering partitions points into k
groups, such that the sum of squares from points to the assigned cluster centers is minimized.
## * No rows dropped by 'na.omit' call
# print our model fit
kmeans_model
## K-means clustering with 3 clusters
##
## Cluster centers:
## Petal_Width Petal_Length
## 1 1.359259 4.292593
## 2 0.246000 1.462000
## 3 2.047826 5.626087
##
## Within Set Sum of Squared Errors = 31.41289
# predict the associated class
predicted <- sdf_predict(kmeans_model, iris_tbl) %>%
collect
table(predicted$Species, predicted$prediction)
##
## 0 1 2
## setosa 0 50 0
## versicolor 48 0 2
## virginica 6 0 44
# plot cluster membership
sdf_predict(kmeans_model) %>%
collect() %>%
ggplot(aes(Petal_Length, Petal_Width)) +
geom_point(aes(Petal_Width, Petal_Length, col = factor(prediction + 1)),
size = 2, alpha = 0.5) +
geom_point(data = kmeans_model$centers, aes(Petal_Width, Petal_Length),
col = scales::muted(c("red", "green", "blue")),
pch = 'x', size = 12) +
scale_color_discrete(name = "Predicted Cluster",
labels = paste("Cluster", 1:3)) +
labs(
x = "Petal Length",
y = "Petal Width",
title = "K-Means Clustering",
subtitle = "Use Spark.ML to predict cluster membership with the iris dataset."
)
Use Spark’s linear regression to model the linear relationship between a response variable and one or more explanatory variables.
lm_model <- iris_tbl %>%
select(Petal_Width, Petal_Length) %>%
ml_linear_regression(Petal_Length ~ Petal_Width)
## * No rows dropped by 'na.omit' call
iris_tbl %>%
select(Petal_Width, Petal_Length) %>%
collect %>%
ggplot(aes(Petal_Length, Petal_Width)) +
geom_point(aes(Petal_Width, Petal_Length), size = 2, alpha = 0.5) +
geom_abline(aes(slope = coef(lm_model)[["Petal_Width"]],
intercept = coef(lm_model)[["(Intercept)"]]),
color = "red") +
labs(
x = "Petal Width",
y = "Petal Length",
title = "Linear Regression: Petal Length ~ Petal Width",
subtitle = "Use Spark.ML linear regression to predict petal length as a function of petal width."
)
Use Spark’s logistic regression to perform logistic regression, modeling a binary outcome as a function of one or more explanatory variables.
# Prepare beaver dataset
beaver <- beaver2
beaver$activ <- factor(beaver$activ, labels = c("Non-Active", "Active"))
copy_to(sc, beaver, "beaver")
## # Source: table<beaver> [?? x 4]
## # Database: spark_connection
## day time temp activ
## <dbl> <dbl> <dbl> <chr>
## 1 307 930 36.58 Non-Active
## 2 307 940 36.73 Non-Active
## 3 307 950 36.93 Non-Active
## 4 307 1000 37.15 Non-Active
## 5 307 1010 37.23 Non-Active
## 6 307 1020 37.24 Non-Active
## 7 307 1030 37.24 Non-Active
## 8 307 1040 36.90 Non-Active
## 9 307 1050 36.95 Non-Active
## 10 307 1100 36.89 Non-Active
## # ... with 90 more rows
beaver_tbl <- tbl(sc, "beaver")
glm_model <- beaver_tbl %>%
mutate(binary_response = as.numeric(activ == "Active")) %>%
ml_logistic_regression(binary_response ~ temp)
## * No rows dropped by 'na.omit' call
glm_model
## Call: binary_response ~ temp
##
## Coefficients:
## (Intercept) temp
## -550.52331 14.69184
Use Spark’s Principal Components Analysis (PCA) to perform dimensionality reduction. PCA is a statistical method to find a rotation such that the first coordinate has the largest variance possible, and each succeeding coordinate in turn has the largest variance possible.
## * No rows dropped by 'na.omit' call
print(pca_model)
## Explained variance:
##
## PC1 PC2 PC3 PC4
## 0.924618723 0.053066483 0.017102610 0.005212184
##
## Rotation:
## PC1 PC2 PC3 PC4
## Sepal_Length -0.36138659 -0.65658877 0.58202985 0.3154872
## Sepal_Width 0.08452251 -0.73016143 -0.59791083 -0.3197231
## Petal_Length -0.85667061 0.17337266 -0.07623608 -0.4798390
## Petal_Width -0.35828920 0.07548102 -0.54583143 0.7536574
Use Spark’s Random Forest to perform regression or multiclass classification.
rf_model <- iris_tbl %>%
ml_random_forest(Species ~ Petal_Length + Petal_Width, type = "classification")
## * No rows dropped by 'na.omit' call
rf_predict <- sdf_predict(rf_model, iris_tbl) %>%
ft_string_indexer("Species", "Species_idx") %>%
collect
table(rf_predict$Species_idx, rf_predict$prediction)
##
## 0 1 2
## 0 49 1 0
## 1 0 50 0
## 2 0 0 50
Split a Spark DataFrame into training, test datasets.
partitions <- tbl(sc, "iris") %>%
sdf_partition(training = 0.75, test = 0.25, seed = 1099)
fit <- partitions$training %>%
ml_linear_regression(Petal_Length ~ Petal_Width)
## * No rows dropped by 'na.omit' call
estimate_mse <- function(df){
sdf_predict(fit, df) %>%
mutate(resid = Petal_Length - prediction) %>%
summarize(mse = mean(resid ^ 2)) %>%
collect
}
sapply(partitions, estimate_mse)
## $training.mse
## [1] 0.2374596
##
## $test.mse
## [1] 0.1898848
Use ft_string_indexer
and ft_index_to_string
to convert a character column into a numeric column and back again.
ft_string2idx <- iris_tbl %>%
ft_string_indexer("Species", "Species_idx") %>%
ft_index_to_string("Species_idx", "Species_remap") %>%
collect
table(ft_string2idx$Species, ft_string2idx$Species_remap)
##
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 50 0
## virginica 0 0 50
sdf_mutate is provided as a helper function, to allow you to use feature transformers. For example, the previous code snippet could have been written as:
ft_string2idx <- iris_tbl %>%
sdf_mutate(Species_idx = ft_string_indexer(Species)) %>%
sdf_mutate(Species_remap = ft_index_to_string(Species_idx)) %>%
collect
ft_string2idx %>%
select(Species, Species_idx, Species_remap) %>%
distinct
## # A tibble: 3 x 3
## Species Species_idx Species_remap
## <chr> <dbl> <chr>
## 1 setosa 2 setosa
## 2 versicolor 0 versicolor
## 3 virginica 1 virginica
Let’s walk through a simple example to demonstrate the use of Spark’s machine learning algorithms within R. We’ll use ml_linear_regression to fit a linear regression model. Using the built-in mtcars
dataset, we’ll try to predict a car’s fuel consumption (mpg
) based on its weight (wt
), and the number of cylinders the engine contains (cyl
).
First, we will copy the mtcars
dataset into Spark.
mtcars_tbl <- copy_to(sc, mtcars, "mtcars")
Transform the data with Spark SQL, feature transformers, and DataFrame functions.
Then fit a linear model using spark ML. Model MPG as a function of weight and cylinders.
# transform our data set, and then partition into 'training', 'test'
partitions <- mtcars_tbl %>%
filter(hp >= 100) %>%
sdf_mutate(cyl8 = ft_bucketizer(cyl, c(0,8,12))) %>%
sdf_partition(training = 0.5, test = 0.5, seed = 888)
# fit a linear mdoel to the training dataset
fit <- partitions$training %>%
ml_linear_regression(mpg ~ wt + cyl)
## * No rows dropped by 'na.omit' call
# summarize the model
summary(fit)
## Call: ml_linear_regression(., mpg ~ wt + cyl)
##
## Deviance Residuals::
## Min 1Q Median 3Q Max
## -2.0947 -1.2747 -0.1129 1.0876 2.2185
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.79558 2.67240 12.6462 4.92e-07 ***
## wt -1.59625 0.73729 -2.1650 0.05859 .
## cyl -1.58036 0.49670 -3.1817 0.01115 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-Squared: 0.8267
## Root Mean Squared Error: 1.437
The summary()
suggests that our model is a fairly good fit, and that both a cars weight, as well as the number of cylinders in its engine, will be powerful predictors of its average fuel consumption. (The model suggests that, on average, heavier cars consume more fuel.)
Let’s use our Spark model fit to predict the average fuel consumption on our test data set, and compare the predicted response with the true measured fuel consumption. We’ll build a simple ggplot2 plot that will allow us to inspect the quality of our predictions.
# Score the data
pred <- sdf_predict(fit, partitions$test) %>%
collect
# Plot the predicted versus actual mpg
ggplot(pred, aes(x = mpg, y = prediction)) +
geom_abline(lty = "dashed", col = "red") +
geom_point() +
theme(plot.title = element_text(hjust = 0.5)) +
coord_fixed(ratio = 1) +
labs(
x = "Actual Fuel Consumption",
y = "Predicted Fuel Consumption",
title = "Predicted vs. Actual Fuel Consumption"
)
Although simple, our model appears to do a fairly good job of predicting a car’s average fuel consumption.
As you can see, we can easily and effectively combine feature transformers, machine learning algorithms, and Spark DataFrame functions into a complete analysis with Spark and R.


TABLE OF CONTENTS
A. GS1 Related Culture 1. ‘Flavours of Bengal in Durga Puja celebrations at Indirapuram Prantic Cultural Society organising five-day event showcasing Bengali pride and heritage; musical performances, cooking competitions and other activities scheduled’ B. GS2 Related Health and Human Resource 1. Child mortality on the decline, says study C. GS3 Related Infrastructure 1. Japan to fund mass rapid transit systems in Gujarat, Haryana Internal Security 1. Manipur tightens security over possible Rohingya influx Energy Security 1. India third in nuclear power installations: study Science and Technology 1. Drug-resistant bacteria emerge from drug companies’ untreated effluents D. GS4 Related E. Prelims Fact F. UPSC Prelims Practice Questions G. UPSC Mains Practice Questions
Context:
- A Bengal-like setting of Durga puja can be seen in Indirapuram, Ghaziabad, where Prantic Cultural Society is organising a five-day cultural celebration starting Sunday. The event will continue till navami, September 29.
Theme:
- ‘Our Heritage, Our Pride’, with an aim to showcase the Bengali pride and traditions associated with the festival.
Keypoints:
- Prantic Cultural society is a socio-cultural forum consisting of groups of people who are passionate to preserve and promote Bengali art and culture.
- It’s a five-day celebration, including musical performances, cooking competitions and theatre by children.
- No artificial colours will be used in the making of the idol.
Category: HEALTH AND HUMAN RESOURCE
1. Child mortality on the decline, says study
Context:
- India’s child mortality rate per thousand live births has fallen by 62% from 125 per thousand live births in 1990 to 47 per thousand live births in 2015.
- This is slightly less than the 2015 Millennium Development Goal of a 66% reduction.
Keypoints:
- ‘Tamil Nadu, Karnataka and Maharashtra have shown remarkable progress between 2000 and 2015’
- The study shows a significant decline in cause-specific child mortality rates between 2000 and 2015 in the country.
- The faster declines in child mortality after 2005 (average annual decline of 3.4% for neonatal mortality and 5.9% for 1 to 59-month mortality) suggest that the country has avoided about one million more child deaths compared to the rates of progress in 2000–2005.
Premature births
- Deaths due to premature births or low birth weight rose from 12.3 per 1000 live births in 2000 to 14.3 per 1000 live births in 2015.
- The low birth weight in poorer States and rural areas was the prime reason behind the increase in the past.
- The Million Death Study titled “Changes in cause-specific neonatal and 1–59 month child mortality in India from 2000 to 2015: a nationally representative survey’ was published online by the Lancet.
- It surveyed the progress in three States during this period, the study says that “if all States of India had achieved the declines seen in Tamil Nadu, Karnataka, and Maharashtra, nearly all States of India would have met the 2015 Millennium Development Goals.”
1. Japan to fund mass rapid transit systems in Gujarat, Haryana
Context:
- $4.5 billion soft loan from Japan International Cooperation Agency to boost $100 bn Delhi-Mumbai Industrial Corridor project
Keypoints:
- The DMIC spans six States (Uttar Pradesh, Delhi National Capital Region, Haryana, Rajasthan, Gujarat and Maharashtra).
- It uses ‘the 1,500-km-long, high-capacity western Dedicated Railway Freight Corridor (DFC) as the backbone’ and aims to be ‘a global manufacturing and investment destination’.
- A soft loan will be utilised to improve two Mass Rapid Transit Systems (MRTS) — one each in Gujarat and Haryana that will be part of the DMIC.
- The JICA is the Japanese governmental agency in charge of implementation of Japan’s Official Development Assistance (ODA).
- Its prime objective is promoting economic development and welfare in developing countries.
- The interest rate of the loan will be kept ‘very low’ (at 0.1%) and have a ‘long’ repayment period (at 40 years, including a 10-year grace period).
1. Manipur tightens security over possible Rohingya influx
Context:
- Manipur has taken stringent security measures at the border areas of Manipur, sensing the possibility of a mass influx in the backdrop of the recent large-scale violence in Rakhine province of neighbouring Myanmar.
Keypoints:
- Assam-Jiribam is the most vulnerable point with respect to the influx of Bangladeshi immigrant, hence additional forces have been stationed at the border.
- Similarly, the drive against illegal immigrants have been strengthened at the Indo-Myanmar border at Moreh by deputing officers.
- The Indo-Myanmar border at Beihang in Churachandpur district would be developed as second corridor in the wake of the Act East Policy adding that a permanent police station would soon be established at Beihang.
- A total ban on use of agriculture land for commercial purpose, “State government is going to effectively implement the Manipur Conservation of Paddy Land and Wetland Act, 2014 by putting a total ban on sale, purchase and use of agricultural land for commercial purposes.”
- Manipuri government reaffirming to effectively implement the Conservation Act 2014,” to keep the records of paddy lands by conducting proper identification of the status of the existing paddy land through digitised cadastral maps on the basis of satellite images to conserve the wetlands.
1. India third in nuclear power installations: study
Context:
- India is third in the world in the number of nuclear reactors being installed, while China is leading at 20, as per the World Nuclear Industry Status Report 2017
- But share of energy generation stagnates globally, and several countries shut down nuclear reactors in 2017
Keypoints:
- Most nuclear reactor constructions are behind schedule, with delays resulting in increase in project costs and delay in power generation.
Nuclear Power vs Renewable Energy Development.
- Since 1997, worldwide renewable energy has produced four times as many new kilowatt-hours of electricity than nuclear power
- The world no longer needs to build nuclear power plants to avoid climate change and certainly not to save money.
- The global nuclear power generation increased by 1.4% in 2016 due to a 23% increase in China, though the share of nuclear energy in electricity generation stagnated at 10.5%.
- By comparison, globally, wind power output grew by 16% and solar power by 30%. Renewables represented 62% of global power generating capacity additions.
Financial crisis
- Massive losses over its nuclear construction projects, Toshiba filed for bankruptcy of its U.S. subsidiary Westinghouse, the largest nuclear power builder in history. AREVA has accumulated $12.3 billion in losses over the past six years.
- French bailout
- The French government has provided a $5.3 billion bailout and continues its break-up strategy.
- Fukushima nuclear power project in Japan, six years after the disaster began, the report notes how the total official cost estimate for the catastrophe doubled to $200 billion.
Category: SCIENCE AND TECHNOLOGY
1. Drug-resistant bacteria emerge from drug companies’ untreated effluents
Context:
- Discharge of untreated effluents from pharmaceutical companies in and around water bodies is contaminated with antibiotics and antifungal agents leading to the emergence of multidrug-resistant bacteria.
Keypoints:
- The bacteria have been found to produce enzymes such as extended spectrum beta-lactamases and carbapenemases which can protect them from antibiotics such as penicillin, cephalosporins, cephamycins, and carbapenems.
- The other antiinfectives found in the waters included antifungal medicine voriconazole, medications for bacterial infections such as moxifloxacin, linezolid, levofloxacin, clarithromycin, ciprofloxacin, ampicillin, doxycycline, trimethoprim and sulfamethoxazole.
- The samples of bacteria tested contain drug-resistant genes.
Nothing here for Today!!!
Nothing here for Today!!!
F. Practice Questions for UPSC Prelims Exam
Question 1. Consider the following phenomenon:
- Fusion is the splitting of a heavy, unstable nucleus into two lighter nuclei,
- Fission is the process where two light nuclei combine together releasing vast amounts of energy.
Select the correct option
- 1 only
- 2 only
- Both 1 and 2
- None of the above
Question 2. Consider the following statements:
- Bharatnatyam dance is known to be ekaharya
- The style was kept alive by the devadasis (performed music and dance as offerings to the deities, in the temple courtyards).
Select the correct option
- 1 only
- 2 only
- Both 1 and 2
- None of the above
Question 3. Consider the following statements
- Child mortality, also known as under-5 mortality or child death, refers to the death of infants and children under the age of five.
- India has the third highest child mortality rate among SAARC countries with 48 deaths reported per 1,000 live births
- India’s child mortality rates declining aggressively
Select the incorrect option/s
- 1 only
- 2 only
- 3 only
- None of the above
G. UPSC Mains Practice Questions
GS Paper II
- The characteristic secrecy in governance is a threat to the objective of “Minimum government and Maximum governance”. Discuss. What effective measures can be taken to set right the path?
- The reducing child mortality rates indicate the commitment towards the Sustainable Development goals. But the study indicates selective performances by the states like Karnataka, Tamil Nadu and Maharashtra. What measures are required by the Central and State governments to achieve parity in development?
GS Paper III
- It is seen that the untreated effluents are the reason behind drug resistance and also health related effects. What possible amendments need to be incorporated in our laws to tackle the same?Comment.
Also, check previous Daily News Analysis
“Proper Current Affairs preparation is the key to success in the UPSC- Civil Services Examination. We have now launched a comprehensive ‘Current Affairs Webinar’. Limited seats available. Click here to Know More.”
Enroll for India’s Largest All-India Test Series


{% extends 'template.html' %} {% block body %} {% csrf_token %}
Your order number is
Product Info: {{ product.title }}
Product Description: {{ product.description }}
Product Price: ${{product.price}}
Continue to homepage


I'm a 23 year old male with an average build. I'm relatively experienced with psychoactives: I've smoked marijuana more or less every day for the last seven years, eaten mushrooms dozens of times, and have limited to moderate experience with other substances that aren't necessary to list. Most importantly to this report, I've smoked Kanna around ten times, both by itself and with other substances - and I know what too much caffeine feels like.
I'd purchased the mesembrine tincture approximately a year ago from an in-country (i.e. canadian) vendor. They said it was 5mg/mL and I trust them, meaning that my 10mL bottle contained 50mg total. I'm not sure what the solvent was, but I'm guessing water. Directly after purchasing it I had tried about 20% of the bottle mixed with a glass of juice. This produced no noticeable effect, except for the fact that I was now down to 40 mg. This then sat on a cool dark shelf in a basement during the intervening time.
Okay, so, the experience. The day in question was yesterday, a somewhat ordinary monday. My boss was out of town, and the nature of my work means that I dictate my own hours, so I'd woken up late morning and was preparing to head off to work around 1 p.m. I love smoking weed and being stoned, lifted, or otherwise out of sorts at work, although I've been smoking a fair amount recently and have started to notice the negative effects I get with a couple weeks of heavy use i.e. general tiredness, laziness, and apathy. Both of these probably factored into my decision to lay off weed for the day and try the mesembrine extract.
At noon-ish (T+0), I put slightly less than half of the bottle (~15 mg) into a glass of water and downed it on an empty stomach. Within a minute I felt a significant mood lift. I took this to be a placebo coming from the knowledge that I was consuming an intoxicating substance, although in retrospect it may have indeed been chemically caused by the mesembrine. If this is the case, the oral absorption of mesembrine is very quick. This decreased over five minutes to just a 'good mood'. I got in my car and ran a twenty minute errand,
[Erowid Note: Driving while intoxicated or tripping is dangerous and irresponsible because it endangers other people. Don't do it!] belting out Desmond Dekker at the top of my lungs the whole way (justification for driving on this substance: I'd smoked Kanna numerous times and never felt it affect me in any significant way aside from slight anxiety relief).
When I returned, I puttered about for a half hour or so, in the back of my head wondering if and when further effects would surface.  Within a minute of the second dose I definitely felt altered. I sat down in front of the TV. A program on Canadian poetry entirely captivated me; this is not unusual for me, but the degree was definitely increased, and I noticed my appreciation for it had at least some of the 'Whoooah, dude, no way. I'm so high' quality. The feeling is hard to describe, although it felt much like coming up on any new substance: slightly anxious and jittery. After about five minutes (T+1hr) the effects ceased to build.
For the rest of the day, I was left with the anxious and jittery feeling. This is not at all what I expected, as Kanna and its alkaloids are noted as being *anti*-anxiety, and my previous experiences were all in line with the slight mood lift I felt at T+0. The significant mood lift I expected with the increased dosage did not materialize. Instead, I felt alert and anxious, but at the same time found it difficult to concentrate. Imagine someone staring straight ahead, wide eyes and open mouth, nervous movements: that was me. Also notable was that my stomach felt irritated, although without any kind of intestinal upset, and I had a slight headache. This slowly decreased over the afternoon, and by five I felt more or less calm, although a co-worker commented that I looked 'confused' - definitely not an ordinary comment. I considered canceling my evening plans but in the end didn't have to.
So, conclusions.
1. Wikipedia lists Kanna as ~0.5% mesembrine, so 40mg mesembrine is equivalent to 80g of Kanna - an unusually large amount. I should do the math beforehand.
2. The closest analogy that I can think of for this experience is a caffeine overdose - that unnecessary anxious feeling that 'something is wrong', coupled with mild body aches. At low doses mesembrine removes anxiety, whereas at high doses it increases it in a fashion similar to caffeine.
3. Depending on your interpretation, the experience might suggest that mesembrine is absorbed very quickly by the stomach: I noticed distinct alterations immediately after both times I consumed it. However, it's also possible that the first alteration was a placebo, and the second alteration was the coincidental come-on of the first dose, which would suggest that the come-on time for oral mesembrine is ~45min. I don't know how to distinguish which of these (if either) are correct.
4. Small doses of mesembrine are enjoyable. Large ones ain't. I'll remember that.
Experience Reports are the writings and opinions of the individual authors who submit them.
Some of the activities described are dangerous and/or illegal and none are recommended by Erowid.


| a1 | a23 | a4 | a56 | ||
| b2 | b3 | b5 | b6 | ||
| c1 | c2 | c3 | c4 | c5 | c6 |
| Reviewed by: ________________________________ | Acceptance: ________________________________ |
| (Producer name / date) | (Representative name / date) |
| Studio | Animes | Last Cell | |
| Pixar | Last Cell 2 | ||
| Studio Ghibli | Grave of the Fireflies | Spirited Away |


{%extends "nav_account.html"%} {%block title%}Meetings{%endblock%} {%block content%}
Your membership is not active, so you cannot participate in any meetings.
{%else%} {%if meetings%}The following upcoming meetings exist. If you wish to attend a meeting, click the link for it which will generate a key for you (if one is not already generated), and then follow the instructions on that page.
There are currently no upcoming meetings registered.
{%endif%} {%endif%} {%endblock%}


Economy-Cars
SUV-Cars
Luxury-Cars
Luxurious Cars
Brand: Aston Martin
Model: DB11
Year: 2019
Type: Luxury
Brand: Bentley
Model: Bentayga
Year: 2020
Type: Luxury
Brand: BMW
Model: M8 Gran Coupe
Year: 2018
Type: Luxury
Brand: Mercedes
Model: S63 Cabriolet
Year: 2020
Type: Luxury
Brand: Porsche
Model: 911 Cabriolet
Year: 2021
Type: Luxury
Brand: Chevrolet
Model: Corvette 69 Stingray
Year: 2020
Type: Luxury


Password filter hooking with no reboot, reflective load/inject DLL
clymb3r/Misc-Windows-Hacking/HookPasswordChange


About/Join
Kottu is a Sri Lankan blog aggregator. It basically collects a slice of the Sri Lankan blogosphere in one place. The only criteria for joining Kottu is:
- Having a working feed (which most blogs do)
- Being ‘Sri Lankan’, as in based in or covering Sri Lankan experiences
- Being original content (not copy/pasted)
- Observing very basic standards of libel and obscenity
- Being updated in the last two month
To join just send a mail to<EMAIL_ADDRESS>I’ll get around to a contact form someday.
What Is Kottu?
 If you’re unfamiliar, a feed is something a blog has. If you’re unfamiliar with a blog, click some of the links on Kottu. It’s basically a personal webpage updated frequently. Every modern blog has a feed. As a blogger this can be irrelevant, but feeds are a standard that allow you to play with data in a lot of fun ways. Basically, you can have hundreds of blogs that look and behave different, but the feeds are all standard. That means Kottu, as a machine, can read them all.
Kottu is, physically, a rented Rackspace Cloud server running Linux and a basic set
of software. One program it runs is WordPress, which is a blogging software itself.
This makes Kottu a sort of metablog. Instead of a writer, however, Kottu reads hundreds
of feeds and posts summaries. You can click on the links to read the original content
on those sites. It also ‘reads’ the feeds from Flickr to find Sri Lankan photos.
There is a lot of technical wizardry behind this, but that is all managed by a
WordPress plugin called FeedWordpress. Setting up Kottu from scratch (now) takes
about 15 minutes and is pretty easy to replicate. Kottu involves almost no original
code, but rather putting existing pieces together.
Who Is Kottu?
Kottu – the name and spur – came from Mahangu. The data structure was imagined by Indi Samarajiva, but it’s not especially originally. It is actually something called a ‘Planet’, used by geeks for years. Mahangu is in Uni and not much in the Blogging scene, but Indi is still an active blogger at www.indi.ca. Kottu was recoded in August 2011 by Janith Leanage and he's now the main Kottu cooker.
Indi, that is I, am basically the janitor. I add feeds, I keep an eye on the database, and generally try to keep the place clean. I do have opinions – which I write on indi.ca – but they have no bearing on Kottu. I’m not saint, but I am lazy and it is easy for me to simply syndicate everyone that fits the requirements above without thinking too much about it. I’m sorry that there’s still a human element here, but there is. I cannot be assed to be any sort of active figure, I’m just trying to keep the place clean.
Concerning the requirements:
- Valid Feed: Basically if you don’t have an RSS or Atom or RDF feed, Kottu can’t syndicate you. It just technically doesn’t work. You don’t need to know what these things are, they come with any modern blog (as listed below).
- Sri Lankan Content: Kottu isn’t a random aggregator, it’s Sri Lankan. That can mean you’re in Sri Lanka, or abroad talking about it, or whatever. This also, by extension, means that overly technical blogs won’t be syndicated because they’re not exactly great for a general Sri Lankan audience. This is obviously a grey area, but blogs 95% devoted to coding or gaming probably won’t fit.
- Original Content: Basically, copying and pasting from newspapers or online sources doesn’t count as a blog post. It’s fine copying a quote or paragraphs of interesting text/image/video and then commenting on it, but what Kottu is interested in is your commentary.
- Libel and Obscenity: The blogging scene is pretty diverse, so please try to think of the community. The occasional ‘fuck’ or ’shit’ is fine, but basically stuff that would totally gross my mother out is not cool. That is, admitedly, a pretty obscure standard, but what to do. In terms of libel, it’s basically not cool to personally attack people to the point of (seriously) threatening violence or publishing phone numbers, addresses, and revealing anonymous identities. There is a lot of rough and tumble on Kottu, but there is a grey line.
- Being Updated In The Last Two Months: Kottu polls the feeds, and it’s not efficient to poll dead ones. If you haven’t updated in the last two months, you will be dropped – without notification.
ETC
Basically, as the janitor, I am busy and I do and pay for this in my spare time. I can’t always respond to email fast, and I sometimes forget to respond at all. I can’t fix your blog and I can only go so far towards resolving problems. I care about Kottu, but this is not a public or even private service and you simply can’t expect very good ‘customer service’. I have a day job and family and life and this is something I look after when I can, so please try to be understanding.
At some point I hope to hand over the admin of Kottu and go a bit open-source, but that’s kinda how it is right now. Anyways, welcome and happy blogging.
Recommendations
If you want to start blogging I wholeheartedly recommend wordpress.com. It’s fully hosted and sets up in minutes. Another big service is Blogger, which is very popular and owned by Google. If you want total control, wordpress.org is what runs Kottu and it’s pretty simple and powerful. You do need to rent a server for that. Kottu runs on a Rackspace Cloud Server, but that’s not necessarily a recommendation.


Add new Material modules you intend to use under material.module.ts.
material.module.ts
Find the Material Icon List here.
Run the following command to pull the Docker image: docker pull duluca/angular4-material2-starter
docker pull duluca/angular4-material2-starter
Learn about Angular Flex Layout here.


untangled
about us
resources
activities
Untangle your stress and anxiety through various stress-reducing activities.
Untangled
activities
Daily Journal
Daily Reminder
Grateful
Games
Resources
About Us


{% extends "split/base.html" %} {% block content %} {% if experiments %}
No experiments have been started yet. You need to define them in your code and introduce them to your users.
Check out the documentation for more help getting started.
{% endif %} {% endblock %}


About
Hi! My name is Deon Taljaard. This is my personal site, which I mainly use as my knowledge base. I enjoy learning and sharing my knowledge. I hope you find it insightful and interesting.
I’m an application developer based in Centurion, South Africa. To develop software is both my job and my hobby and I love it.
I believe in having a balance in life. Thus, in my personal time, there are a host of things that I keep myself busy with. Amongst other things, I:
- serve as youth mentor / leader at our congregation
- spend time with my fiancée
- hang out with friends
- work on a personal project with my best friend
- try (yes, try) and stay fit in the gym
- play stuff on my acoustic guitar
- lehre mich Deutsch
I never get time to do all of the above in one day, however, I never sit around twirling my fingers and try to make the most of the days on Earth I have been blessed with.
A big thank you to Bronwyn Alsyon Photography for the profile pic and logo design.
Thanks for reading!


Simple Easy
From plug in to powering on the ease of use for Bubble is entirely seamless. Within a few seconds you can be connected to Bubble via WiFi.
Any Device
Laptops, tablets, and phones alike can access any content you wish to be able to distribute. Bubble can stream media content simultaneously to any device that supports modern browers, such as Safari, Chrome, Internet Explorer and Firefox.
Ultra Slim
Bubble is pocket-size portable with dimensions of 0.98"H x 2.5"W x 3.6"L. It has an ultra light yet durable construction. And let's not forget it's entirely Internet free.


Class MD5, for calculating MD5 hashes of strings or byte arrays it is not meant to be fast or secure.
More...
#include <md5.h>
Class MD5, for calculating MD5 hashes of strings or byte arrays it is not meant to be fast or secure.
usage: 1) feed it blocks of uchars with update() 2) finalize() 3) get hexdigest() string or MD5(std::string).hexdigest()
assumes that char is 8 bit and int is 32 bit
◆ size_type
◆ uint1
| typedef unsigned char uint1 |
◆ MD5() [1/2]
◆ MD5() [2/2]
| MD5 |
( |
const std::string & |
text | ) |
|
◆ finalize()
◆ hexdigest()
| std::string hexdigest |
( |
| ) |
const |
◆ init()
◆ update() [1/2]
| void update |
( |
const unsigned char * |
buf, |
|
|
size_type |
length |
|
) |
| |
◆ update() [2/2]
| void update |
( |
const char * |
buf, |
|
|
size_type |
length |
|
) |
| |
◆ operator<<
| std::ostream& operator<< |
( |
std::ostream & |
, |
|
|
MD5 & |
md5 |
|
) |
| |
|
friend |
◆ digest
The documentation for this class was generated from the following file:
- /home/ricardo/workspace/curro/eProsima/desarrollo/proyectos/fastrtps/include/fastrtps/utils/md5.h


{% extends "_base.html" %} {% block content %} {% if current_user.is_authenticated %}
Name: {{users.user_firstname}} {{users.user_lastname}}
Email: {{users.email_user}}
Click here to logout: Logout
If you would like to delete the information and account, please click: Delete Account
{% endif %} {% endblock %}


{{ comment }} This page is for showing the business associates that an entity has an agreement with. If an associate is also a covered entity, it has a link to more information regarding it. Each agreement with a busines associate has information that the entity and business associate share so the business associate can fulfill the agreement. If information access is restricted to the accessor, the information shared cell can contain the reason the information is being restricted. {{ endcomment }} {% extends "base.html" %} {% block title %} Business Associates {% endblock %} {% block top %}
| Entity {% include "invisible.html" %} | Information Shared {% include "invisible.html" %} | Purpose {% include "invisible.html" %} |
| {% if association.BusinessAssociate.Covered %} {{ association.BusinessAssociate.Name }} {% else %} {{ association.BusinessAssociate.Name }} {% endif %} | {% with data=association.SharedInformation %} … {% include "informationModal.html" %} {% endwith %} | {{ association.Purpose }} |